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CHAPTER 1 


INTRODUCTION 

Most of the commercial alloys used for several 
structural components have two-phase or multi-phase structures. 
Examples of such alloys are plain carbon and alloy steels, nickel 
base super alloys, brasses, titanium alloys and aluminium alloys. 
Among steels the specific examples are spheroidized steels 
( f err 1 te-cement i te ) , dual phase steels (ferrite-martensite) and 
duplex stainless steels (austenite-ferrite). Alpha+beta titanium 
alloys such as Ti-6A1-4V are examples of alloys having B.C.C. ft 
phase and H.C.P. ot phase. Many advanced materials which are 
currently being developed, like i ntermetal 1 i cs and metal matrix 
composites, also have two-phase/mul t i -phase structures. The 
presence of more than one phase often enables the development of 
structures which are efficient in load carrying applications and 
their properties can be varied over a wide range [ID. 

1.1 Mechanical Properties of Two-phase Alloy 

Mechanical properties like y i e ld-strengthC 1 -53 , 
hardness[13, duct i 1 i ty[ 53 , creep[43 , fatigueC53, fracture 
modeC53, and the coefficients of constitutive equation describing 
plastic flow£5,63 in a given alloy have been found to vary 
significantly with the distribution and volume fraction of phases 
present. This dependence is over and above their variation with 
grain size, sub-grain size. texture, coherency strain between 
different phases and interface phase, if present. 
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Fig.(l.l) shows the typical variation of yield strength 
versus the volume fraction of second phase in two phase 
alloysCl]. Ankem and Margolin[23 have given the following 
expression for the yield strength of a two phase alloy. 


o = 

Y 


va 


+ °vft V /3 + l aft 


. . . <1 .1 ) 


where V and are volume fractions of a and ft phase, a and 

aft va 

G V 

are their bulk yield strengths and the term 1^ is due to 

the interaction between the phases. This expression shows that 

the law of mixture for the yield strength need not be followed. 

The yield strength of the alloy not only depends on the volume 

fractions of two phases, V and but also on the interaction 

oi fy 

Y 

parameter 1^ which is likely to depend on the morphology of two 
phases . 

Similarly, hardness for two-phase alloys has been 
shown to vary asCID, 


H = H^d - CV^> + H^tCV^) ...<1.2) 

where H and are hardnesses of a and ft phases, is the 

aft ' ft 

volume fraction of /9-phase in the alloy and C is termed as 
cont i gui ty . 

Fig. <1.2) shows true-stress vs true-strain curves for 
four a+ft Ti alloys with different volume fraction of a-phaseC53. 
Significant rise in stress-value at a given value of strain (much 
above yield strain) with volume fraction of a phase can be 
observed. Sastry et alC63 observed that flow stress as a function 
of strain-rate, significantly varies with volume fraction of 


a-phase . 



0.8 

VOLUME FRACTION OF /3 


Fig.l 


Fig. 1.2 


. 1 Typical curves showing various forms of the 
variation of yield strength versus volume 
fraction ft of two-phase alloys 



True stress - true strain curves for four 
(a+ft) Ti alloys having fine a-phase particle 
and different volume fraction of a-phase 
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A sign if icant amount of work has been done on the 
variation of piro*>e rt i es of ti tanium al loys with phase morphology. 
Tablefl.l) shouws the variation in tensile properties of Ti-6A1-4V 
alloy with different thermo-mechanical treatments i.e. with 
d if fe rent sizes san d morphology of c*-plates. It can toe seen that 
both yield strength and ductility can be varied significantly. 
Similarly, Tab le <1 ,2) shows that fatigue properties also vary 
subst ant i al ly uwi ~th ex-plate morphology. Creep strength is also 
sensi ti ve to a — p la te morphology in Ti-alloys and is the basis of 
evolving modif ie«d heat treatments in several high— temperature 
titanium alloys, Table <1.5) shows the variation of toughness in 
T i-6A 1-4V allo^ uwi th ot-plate morphology. M i c rostructun 1 features 
thus play a vital role in mechanical properties of 
two-phase/mult :i-*phase alloys. 

1.2 Objectives of Designing a Deformation Processing Schedule 
for Two — pbiase Alloys 

One of the most important challenges -Faced by the 
design engineer's and manufacturers of various commercial alloys 
is obtaining thie desired microstructure in the as processed 
material so that the desired property levels in them can be 
obtained. This is conventionally done by utilizing the 
heat- tr eatment principles. However. the recent increase in 
emphasis on cdarwage -tol erant designC33 has led to the realization 
that mi crostructi_ir e must also be controlled on a scale which is 
not s igni fi canffclv affected by heat treatment . An effect ive way of 
manipulating micro structure on such a scale is through combined 
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TABLE 1 . 1 


Tensile Properties of Three a+/? Ti-Alloy C Ti -6A1 -AVD [ 3] 


SI. 

No. 

Morphology 
of a-phase 

Condition 

Yield 

stren- 

gth 

Tensile 

strength 

Elonga- 

tion 

Reduc— 
tio in 

area 




CMPa) 

CMPa) 

CX> 

OO 


equiaxed 

a+ft forge + 
recrystal 1 i - 
zation anneal 

711 

876 

12.4 

56 

equ iaxed 
(coarse) 

a+ft forge + 
mill anneal 
(mi nm. values) 

828 

897 

10.0 

25 

equiaxed 
+ plate 

1 ike 

a+ft forge + 

STA (age 4h 

at 594°C) 

876 

958 

15.2 

54 

equiaxed 
+ plate 

1 ike 

a+ft forge + 
STOA (age 24h 

at 594°C) 

904 

975 

15.5 

47 

plate like 

ft forge + AC 

CBA> (+705°C/ 
2h/AC> 

775 

856 

11.2 

25 

plate like 

ft forge + WQ 

CBQ> (+705°C/ 
2h/AC) 

865 

952 

5.9 

6 

equiaxed 

a+ft forge DA 
(870°C/2h/AC 
+ Mi° C/ 2h/AC) 

856 

911 

15.5 

47 
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TABLE 1 . 2 


Fatigue and Tensile Data for Various Microstructural 



Conditions 

; of a+/3 Ti-Alloy 

CTi-6Al 

-4VM31 


Si. 

No. 

Condition and 
a— phase morpho- 
logy 

Yield 

stren- 

gth 

CMPa) 

Tensile 

strength 

CMPa) 

El onga- 
tion 

C» 

Stress 

C smooth!) 

at 10 7 
cycles, 

CMPa) 

Stress 
C notc- 
hed) at 

1 0 7 cy- 

cles, 

CMPa) 

1 

10* equiaxed 
primary a + 
anneal 

971 

1068 

14 

537 

214 

2 

40* equiaxed 
primary a + 
anneal 

930 

1013 

15 

579 

225 

3 

10* equiaxed 
primary c*. + 

STOA 

978 

1061 

15 

489 

220 

4 

30* elongated 
primary o + 
anneal 

923 

1020 

13 

620 

227 

5 

ft forge + anneal 
(plate like) 

882 

992 

11 

565 

220 

6 

ft forge + STOA 
(plate like) 

978 

1075 

10 

586 

220 

7 

10* equiaxed 
primary a + 
anneal 

882 

985 

13 

620 

214 
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TABLE 1 . 3 


Fracture - Toughness Variations in a+/9 Ti-Alloy CT1-6A1-4VD 

tor Plate Shaped Specimen! 31 


SI. 

Mor phhol ogy 

Condition 

Yield 

Tensile 

Elonga- 

K CK 3 
xc a 

No. 

of a -phase 


stren- 

strength 

tion 




gth 

CMPa> 

C'HPaD 

CJO 

MPa jin 


i 

equiaxed 

cx+ft roll + 
mill anneal 

1096 

1171 

14 

32 

2 

equ i axed 

a+fi roll + 
recrystal 1 i - 
zation (RA> 

1034 

1144 

13 

51 

3 

equ i axed 
plate like 

a+(3 roll + 
duplex anne- 

882 

971 

14 

59 



al(954°C/l . 5h/ 
AC+760°C/lh/AC) 




4 

plate like 

BA (a+/3 roll 

875 

951 

11 

87 



+1 038°C/20mi n/ 






AC +732°C/2h/ 
AC> 





5 

equiaxed 

a+/3 roll + 

785 

882 

— 

94 


recrystal 1 i - 
zation anneal 
(RA> 
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thermal and mechanical working. Such a processing schedule allows 
a broader range of microstructures which can be obtained than 
that is possible through thermal treatment alone. Such 
manipulations are generally known as thermo mechanical 
process ing(TMP) . 

Objectives of designing a deformation processing 
schedule for two-phase alloys therefore are : 

(a) transformation of initial microstructure from an 
initial to a final microstructure 
<b) assessing load requirements so that a given 
deformation processing schedule can be carried out 
with a given set of equipments 
<c> assuring that the work piece material is workable 
under the given stress - strain - strain rate 
temperature domain. 

The present study is not concerned with the third objective. 

1*2.1 Transformation of the initial microstructure to a 

final one 

Mechanical working has been historically used 
as the primary means of changing the size and shape of materials 
while transforming the cast structure of an ingot into what is 
generally refered to as a wrought product. The two principal 
methods of working are forging and rolling along with other means 
such as extrusion and drawing. But at present it is also done for 
obtaining desired microstructure. Mi crostructural features of 
importance for property control includeC33 : 
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(a) coherency and distribution of strengthening 
prec i p i tates 

<b> degree of recrystallization 

(c> grain size and shape 

(d) crystallographic texture 

<e> size, shape and distribution of dispersed second 
phase <e.g. i ntermetal 1 i cs > . 

Examples of incorporating thermo-mechanical working in 
the manufacturing schedule are controlled rolling of steels, 
thermo-mechanical processing of A1 alloys and control of 
morphology of a-plates in a+/3 Ti alloys. During mechanical 
working of Ti-alloys in the a+(3 range, a-plates elongate with 
deformation. Finally they break into smaller plates which get 
recrystallized during subsequent heat-treatment into equiaxed 
structures. Fig. <1.5) shows the effect of strain due to forging 
and subsequent annealing on the length distribution of 
a-phase[73 . 

Roberts has reviewed the mi crostructural changes during 
hot deformation of coarse two-phase materials like a+ft brass and 
a+r stainless steel[83. By coarse two-phase material we mean that 
both the phases have same grain-size. There have been very few 
systematic studies aimed at elucidating the mi crostructural 
evolution in association with hot deformation of two-phase 
alloys. In these studies it has been found that one 
microconstituent undergoes restoration through dynamic recovery 
plus recrystallization and the other softens via dynamic recovery 
only. It has been observed that orientation of plate like phase 
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Fig. 1.3 The effect of strain by forging and subsequent 
annealing on the length distribution of the 
ot-phase , in material with initial thin 
lenticular structure 
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significantly changes with time during deformation processing. 


1.2.2 Assessment of the load requirement 

For designing a high temperature deformation 
processing schedule under given conditions it is very important 
to know the loads, size reduction and overall rate of deformation 
imparted to a two- phase alloy which will lead to the desired 
microstructure. Many theoretical studies using analyt ical C93 and 
computat ional [10] procedures have been made to meet these 
requirements in rolling, forging, extrusion and other metal 
working operations. These procedures, however, treat the material 
only at macro-level and mi crostructural variables have not been 
accounted for. However, a few efforts have been made in this 
direction in the recent past. 

GurlandCl] has reviewed the application of the law of 
mixtures to the plastic deformation behaviour of two-phase 
alloys. He has shown that the law of mixture which states that. 


a 

c 


O' V + o' V 
ot a ft ft 


...(1.5a) 


e - e V + 
c ot a ft ft 


...<1.5b) 


where V is volume fraction, and a and e are the average values of 
the in-situ stresses and strains in the composite and in each of 
the two-phases respectively, is applicable for the plastic 
deformation of the two-phase alloys in the small strain range. 


Sastry et al[6] have discovered that law of mixture is 


not applicable .for the case of high 
def ormat i on in a+ft Ti alloys. 


temperature 


plast i c 
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Ankem and Margolin[53 by using NASTRAN computer program 
and stress-strain curves for cn and ft phases of T i -Mn alloys have 
calculated reasonably well the effect of particle size, matrix 
and volume fraction on stress-strain relations of the two-phase 
a+ft Ti-alloys. They found that for a given volume fraction of 
a-phase the calculated stress-strain curve was higher for a finer 
particle size than for a coarse particle size and this behaviour 
was seen for all the alloys with different volume fraction of 
a-phase. The calculated stress-strain curves for four a+ft alloys 
with different volume percent ofa-phase were compared with their 
corresponding experimental curves and in general good agreement 
was found. Discrepancies were attributed to lack of consideration 
of particle morphology. 

Dragone and Nix have studied the creep behaviour of 
fiber reinforced metals using FEME11D. They have modelled a 
fiber reinforced metals shown in Fig. (1.4a) whose creep 
deformation rate will depend on many factors like fibre volume 
fraction (V f ) , fibre aspect ratio (L/D) , degree of fibre overlap, 
fibre/matrix interface properties (a. . ), relative fibre 

interface 

spacing (D/L) as simple unit cell shown in Fig. (1.4b). Their 
results are shown in Fig. (1.5). A study of this type can be 
useful for understanding the effect of morphology of individual 
phases on the overall high temperature plastic deformation of 
two-phase alloys. This is very important from processing point of 


view. 
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Fig.l.^a Variation in the calculated composite steady 
state creep rate with fiber aspect ratio at 
constant volume fraction for two different 
cases: <1> fiber aspect ratio varied while 
keeping unit cell aspect ratio <l/d> equal to 
fiber aspect ratio <L/D> and (2) fiber aspect 
ratio varied while keeping relative fiber 
spacing (d/L> constant 



.01 .1 l 10 

Normalized Inter-Reinforcement Distance d/L 

Fig.l.5b Variation in the calculated composite steady state 
creep rate with relative reinforcement spacing for 
various unit cell geometries. Both fiber and plate 
geometries are included on this plot 
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1.3 Limitations of Existing Analytical Procedures 

A review of the literatures on deformation processing 
of two-phase material, therefore, indicates the following 

(a) most of the computational work associated with the 
load requirement during plastic deformation of 
two-phase structures has been carried out assuming 
material to behave like e las to-p last i c . This work, 
though useful, is strictly applicable for low 
temperature deformation 

(b) very little work has been done in connection with 
heterogeneities in stress, strain and strain-rate 

<c> effect of shape of the second phase particles on 
the overall deformation behaviour has not been 
accounted for 

(d) most investigators have analysed the overall 
deformation behaviour of material assuming that 
the second phase plates/fibers are aligned in a 
particular direction. This is applicable for 
aligned fiber composites but not for two-phase 
al loys 

(e> in principle, overall deformation properties 
assigned should be invariant to e, e fields. This 
check has not been carried out by most of the 
earlier investigators. 



CHAPTER 2 


FORMULATION OF THEORETICAL MODEL FOR THE HIGH TEMPERATURE 
PLASTIC DEFORMATION OF A TWO PHASE STRUCTURE 

2.1 Introduction 

As stated in the previous chapter, the response of a 
two-phase structure to the stress state imposed on it during its 
deformation processing depends on the plastic deformation 
behaviour of the aggregate which in turn is influenced by : 

<i> physical and mechanical properties of constituent 
phases 

<ii> their volume fractions, and 

( i i i )morphologi cal features i.e. the shape, size and 
distribution of individual phases. 

Due to the large number of variables which are usually needed 
for obtaining an explicit constitutive relationship for a given 
structure the plastic deformation behaviour of a two-phase 
aggregate is known to be extremely complex. ZaouiC12] has pointed 
out two possible approaches which in principle can be adopted to 
model the plastic behaviour of a given material. These are : 

(a> the first one consists of phenomenological 
approach using internal variables whose physical 
meaning is not to be elucidated. These 
parameters are only expected to allow after 
adequate identification from standard experiments, 
a fairly good predicting ability i.e. fairly good 
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aggrement with further experimental results 
(b) the second one proceeds from a deductive approach. 
It starts from some knowledge of physical 
phenomena involved in the macroscopic response and 
works out a mechanical description of related 
variables and laws of evolution. Then by dealing 
with available information (or adequate 
assumptions) on the initial microstructure of 
concerned material, it has to deduce the 
corresponding macroscopic behaviour through some 
homogenization procedure. 

It is the second approach which has been adopted here. Clearly 
modelling of the deformation behaviour of a two-phase aggregate 
requires an understanding of deformation behaviour of individual 
constituents. The approach for the high temperature deformation 
of two-phase structure which has been adopted in the present 
analysis has been described in this chapter ... Thi s i s pr.ecee.ded by 
a general discussion on the fundamental considerations of the 
plasticity theory. 

2.2 Plastic Deformation 

2.2.1 Stress, strain and strain-rate tensors 

The quantitave description of plastic 
deformation (as a matter of fact any type of deformation) is done 
by relating the stress tensor to strain and strain rate tensors. 
So, before reviewing the phenomenon of plastic deformation, 
stress, strain and strain-rate tensors will be defined. 
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State of stress in three-dimensions is described by a 


tensor of second rank, i.e. it consists of nine components, six 
of which are independent. Thus, 



where a = a 
12 21 


cr 

13 


a and 

31 


a 

23 


a 

32 


Like all other symmetric tensor quantities a., can be split into 

vj 

two components namely (i) hydrostatic and (ii) deviatoric. 
Accordingly, 

a . . = a. . + i a . . 6. . ...(2.2) 

Vj VJ > IJ 

where o'.. is the deviatoric component and 6.. is the Kronecker 

VJ Vj 

delta defined as. 


1 0 0 


6 . . = 0 1 0 
vj 

0 0 1 


1 i = j 

. . . (2.3> 

0 i * j 


It is clear from the above definition that the stress state at a 
point changes the orientation of coordinate axes. Invariant 
coefficients of stress tensor i.e. I , I , I are obtained in the 

12 3 

usual manner. They are as follows. 
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and 

S imi larly f 

£ 

i i 

where s 

xx 

£ 

XV 

£ 

xz 

£ 

YZ 

here 6 , < 

x 

X, Y f and 
with time. 


I = a + a + a 

1 XX YV ZZ 


. . . (2.4a) 


I = 


a a + a a + a a - a 

XX YY YY ZZ XX ZZ XY 


2 2 

o' - a 
XZ YZ 


(2.4b) 


I=o o o +2o o o - o o 

a XX YY ZZ XY YZ ZX XX xz 


o o 

ZZ XY 


. . (2.4c> 


strain tensor at a point is defined as 
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... (2.5) 
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FT * £ S5Z ~ 3z~ 


, 66 66 
1 / X ^ Y . 

£ YX " 7 ( W~ + Wy~ } 


= £ 


ZX 


, 66 66 

1 / * . * 

2 61 6X 


Z Y 


. d<5 d<5 

t ( w + wr~ } 


5 and 6 are components of displacement vector along 

Y Z 

Z axes. Strain-rate is the rate of change of strain 
Therefore, strain-rate tensor is defined as. 
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here U, V and W are velocity components along X, Y and Z 
di rect i ons . 

Strain and strain tensors also like stress tensor can 
be split into two components as hydrostatic and deviatoric. 




t 6. . + e . . 
> u. vj v 4 


. . . (Z.7) 


where e. . is deviatoric component of strain tensor and 

»• j 


T £. . O. . + 6 . . 

> w 14 v 4 


. ( 2 . 8 > 


where « is deviatoric component of strain-rate tensor. 


Invariant coefficients of strain and strain-rate 


tensors are defined in the way similar to that followed for 
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stress-tensor . Effective stress (<?•>, effective strain is) and 
effective strain-rate is) are related to second invariant of 
their respective tensors and are expressed as. 


0 = \ 


* = \ 


and 
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» • • » 2 . - • • v 2 . . • • x 2 . -.,2 . .2 , .2 \ 

<£T -e> + <£ -«) + <« -c)+ 6(« + £ + S ) 

XX YY YY ZZ ZZ XX XY YZ XZ 


. . . (2.9c) 


2.2.2 Incompressibility 

If it is assumed that during the course of plastic 
deformation of any material no cavitation or internal cracks are 
formed, the density of the material remains unaltered. Under such 
conditions the constancy of material volume is maintained. 
Therefore, the sum of all the diagonal components of strain 
tensor is. . ) , which represents volumetric change, is zero. Thus 


= s 


n 


+ s 


22 


+ S 


aa 


= 0 


. . . ( 2 . 10 ) 


This along with Eqs . 2 . 7 and 2.8 lead to 


f 


£. . = £ . . 
V) vj 


. . . <2.11a> 
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and . . . ( 2 . 1 lb> 

This has an important consequence in theories of plasticity which 
so far have been attempting to relate the deviatoric component of 
strain and strain-rate to deviatoric component of stress. 

2.2.3 General constitutive equations 

Plastic deformation in general can be 
expressed by constitutive equations of the type, 

a - f (£ , % , M> ... (2. 12) 

where M represents the microstructural state. 

Most materials below the recrystallization temperature, 
i.e. in the cold working range, are not affected by moderate 
strain rates. Approximate stress-strain relationships for a 
limited region of strain can often be given by exponential 
equation of the form, 

a = K(c> n ... <2.13) 

where K and n are constants. 

Apart from the above relationship other constitutive 
relationships have also been suggested. These are summarized in 
Table 2.1. 

It can be assumed that at high temperatures dynamic 
softening mechanisms occur concurrently and hence no 
work-hardening occurs in the material being deformed. The most 
commonly used expression for high temperature deformation which 
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is the case in the present study is, 

* 

a = KU> m ... (2.14) 

where rr/is the strain-rate sensitivity. The above constitutive 
relationship suggests that the material behaves like a Newtonian 


fluid when 

/ 

m = 

1 . For values 

of rr/s* 1 , the 

material can 

be 

assumed to 

be 

undergoing a 

v i scoplast i c 

flow. Viscoplastic 

treatment of 

plastic deformation behaviour 

of a ma t e r i a 1 

i s 


discussed in sections 2.5 and 2.4.2. 

TABLE 2. 1 

Different Constitutive Relationships for Plastic 
Deformation in the Cold Working Range 


SI. No. 

Author 

Rel at i onshi p 

1 

Unknown 

a = K<e) n 

K and n are constants 

2 

Ludwi k 

cr = a + b<£ )° 

a, b, c are constants 

5 

Voce 

a - a + Cb-a] C 1 -exp( -Ce ) 3 

a, b. c are constants 

4 

Swift 

a = CEa+fi3 n 


2.3 Viscoplasticity 

It is well known that under tensile or compressive load 
a ductile material deforms elastically upto a certain strain 



beyond which it starts yielding. This limiting strain is known as 
elastic limit and upto the elastic limit the stress and strain in 
most of the materials are related by Hooke's law. Beyond the 
elastic limit when the material yields plastically, the stress 
(i.e. flow stress) becomes a function of strain and strain rate. 
Some of the properties of the solid undergoing plastic 
deformation are similar to that of a Non-Newtonian fluid i.e. the 
relationship between the shear stress and shear strain rate is 
non-linear. In fact, it is possible to define a psuedo-v i scos i ty 
for a plastically yielding material which is the ratio between 
the various shear stress and the shear strain rate components. In 
the limit of a perfectly elastic or rigid solid the viscosity 
assumes an infinitely large value. 

2.4 Viscoplastic Deformation Model for a Two Phase Composite 
2.4.1 Definition of the problem 

As mentioned earlier the plastic behaviour of 
two phase aggregate is a complex one. One particular expression 
of, this complexity, strictly reflecting that of nature, consists 
in the large number of variables which are usually needed in 
order to yield an explicit constitutive law which could suffer 
any extensive comparison with experimental data for a given 
mater iaICIZJ . Such a difficulty may be overcome by some 
homogenization procedure. By homogenization we mean finding out 
deformation properties of a homogeneous material whose 
deformation behaviour is same as the overall macroscopic 
deformation behaviour of the two phase aggregate. It can be done 



25 


in either of the two ways : 

(a> microstructure is completely described and the solution 
may be an exact one 

<b> microstructure is only partially and statistically' 
def i ned . 

First case is very difficult and impractical. A limited amount 
of work based on the second approach has been carried out so 
far[12] . In principle it is possible to work out many 

homogenization schemes. 

In the microstructure of many alloys like Ti-alloys as 
well as in aligned fibre composites, plate like second phase 
structures are observed as shown in Fig. (2.1). The structure 
shown in Fig. (2.1) is very similar to that in Fig. <2. 2). The 
deformation behaviour of the two phase aggregate shown in 
Fig.<2.2> can be very useful to pradict the overall macroscopic 
deformation behaviour of two phase aggregate which have plate 
like structure. As a first step, it was decided to study the 
deformation of one of the element of the aggregate. So, in the 
present investigation it was decided that a thorough study of a 
geometrically simple two phase composite, as shown in Fig. (2.5) 
can enlighten many features of two phase deformation. 
Investigation consists of modelling the deformation of a 
rectangular plate-like viscoplastic second phase particle 
embedded in a first phase viscoplastic matrix, within a square 


doma i n . 




Fig. 2.1 Schematic diagram of 
two phase (ct+p) micro- 
structure with plate - 
like oL - phase. 
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Fig. 2. 2 (a-fp) microstructure 
(simplified) 



Fig. 2.3 A unit cell taken from Fig. 2.2 
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2.4.2 Constitutive relation for visco-plastic flow 

Psuedo-V i scos i ty of metals at high 
temperature is dependent on local strain rates. The constitutive 
relation is 

o'.' . = 2 fu s'. . . . .<2. 15) 

»■ j v j 

« 

where e. .is the deviatoric part of strain-rate tensor. The 
1 j 

strain-rate tensor is expressed in terms of velocity gradients in 
the material . 

The viscosity fj for a strain-rate sensitive material is 
taken to be a function of uniaxial yield stress of the material 
<<y > and the effective strain-rate < k ). 

Y 


P = f (O' , k ) 

Y 


1/n 


Cf + [-pr— 3 

Y TW 


•i? * 


. . . <2. 16) 


where y and n are the physical constants which define the 
visco-plastic characteristics of the material. Values of y and n 
usually lie between 1 and 2. 

Flow stress or effective yield stress a of the material 
is given by 


a 


* «- £ ,4/n 

"v + l Tfr 3 


. . . <2.17> 


By rearranging the terms in Eqs.2.16 and 2.17, 
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Or 


a 


■fi P ^ 



But since 


a - 


K e 


Therefore , 




K 


m-± 

r 


...(2.18) 


. . . <2. 19> 


2*4.3 Velocity fields in a deforming visco-plastic 

material 

Constitutive relation for visco-plastic 
materials is given by Eqs.2.15 and 2.19. 

From mass balance (or incompressibility condition) 


or 


-> -> 

V .V = 0 




. . . ( 2 . 20 ) 


for two-dimensional deformation. 
From stress balance. 


-> -> -> -> 
p(V .V .V > = V *«r + f 


where 'f* is the body force and 


<y = <y ii+o- ij+o" j i + <y jj 

= XX XY YX YY 


. . . ( 2 . 21 ) 


. . . < 2 . 22 ) 
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The stress balance equation can be thus written in the form 


9U 


dU 


X direction : p(U + V ^ ) 


= {!>? <0 'xx > + 57 ( 


a )} 

XY J 


. . (2.23a) 


Y direction : p(U + V ) = -[|^<^ xy > + 37 ( ^ YY >} 


. . . (2.23b) 


The constitutive relations for 2-dimensional deformation can be 
written as follows. 


xx 


■p + 2 p 


0U 

W 


XY 


.au . av. 
‘ "<57 + 5x> 


and a - -p + 2 y. 


YY 


£V 

dY 


. . . ( 2 . 24a> 

. . . (2.24b) 

. . . (2.24c) 


O'. . 

where p is pressure or hydrostatic stress (or — y— ) 

Substituting the constitutive relations in Eqs.2.15 and 2.19 we 
get , 


p(U 


p(U 


du 

dU 

w + 

v 37 


u dV 

m + 

v 37 


+ 2p|^) +|y [^<§7 + §*>]}• • • ( 2 . 25 a) 

{§7 <_p + ZfJ W y + fx + • • <2. 25b) 


The above two equations along with Eq.(2.20) form the set of 
governing equations describing the visco-plastic model. 
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2.4.4 Boundary conditions 

Once the governing equations required to 
model the deformation of individual phases have been laid down, 
one needs the boundary conditions also. Boundary conditions that 
are required to solve the governing equations can be classified 
into three categories : 


(a> essential boundary condition i.e., specifying U.P.V 
along the boundaries 

(b) natural boundary condition i.e., specifying ^ 

dV 

and along the boundaries 

<c> traction boundary condition i.e., specifying t and t 

x y 

along the boundaries. 

For the present investigation any of the above type of 
boundary condition can be used. For the sake of simplicity 


essential boundary conditions have been used. 


2. 4. 5 Variable parameters 

The next aspect of the model is the 

characterization of shape, size and relative orientation of 

plate-like second phase particles. Shape of the second phase 

plate-like particles can be conveniently described by L/D ratio 

where L is the length and D is the width. Further L x D gives the 

size of these particles (i.e. volume fraction). 

While simulating this model these three parameters have 

K 

been varied. The ratio where K and K stands for the values 

l\ 12 

1 

of constant K as defined in Eq.<2.14) for the matrix and the 
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var i ed 


plate-like second phase respectively, also have been 
quite extensively. Boundary condition is also one of the variable 


parameters that has been used in the 

present 

i nves t i gat i on . 

The 

value 

of 

m’ has been kept constant for both 

the phases 

and 

set 

equal 

to 

0.05. The range over which 

these 

parameters 

vary 

i s 

given 

i n 

Table 2.2. 






TABLE 2. 2 

Range of Variation of Model Parameters 


SI. No. 

Parameter 

Range 

1 

volume fraction 

( V 

0.08 

- 0.44 

2 

L/D 

1.0 

- 4.0 

3 

K /K 

2 1 

2.0 

- 5.0 

4 

plate orientation 

(m) 

0.0 

- 7.0 


2. 5 Possible Approaches to Solve the Governing Equations 

Analytical solutions for differential equations in 
engineering problems are available for very rare situations. 
Apart from that in the present investigation two phases with 
different material properties are there. It was decided to solve 
the governing equations with given boundary conditions 
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numerically which can be made possible with the help of computer. 
Out of various computational techniques for solving differential 
equations finite element method is one of the most powerful 
technique available and has been utilized in the present study. 
Finite element analysis of the present model is discussed in 


chapter 5. 



CHAPTER 3 


FINITE ELEMENT ANALYSIS 

3.1 Introduction 

As stated in chapter 2 to model the deformation of a 
two-phase alloy at high temperature, the differential Eqs.2.20, 
2.25a and 2.25b, which governs the flow of viscoplastic material 
have to be solved. Finite element technique employing Galerkin's 
weighted residual method used to meet this objective. In this 
chapter first one of the numerical techniques used in our 
investigation, weighted residual method has been discussed. Then 
finite element procedure has been outlined and its application 
for solving the governing differential equations has been 
outlined. At the end of this chapter the essential features of 
the programming have been presented. 

3.2 Weighted Residual Method 

A differential equation should be satisfied at every 

point in the solution domain. However the numerical solution for 

the problems may not satisfy the governing equations exactly. It 

may satisfy the equations approximately having a small non-zero 

residue at most of the regions in the solution domain. 

* 

Let <f> be the exact solution and 4> be an approximate 
solution. Substituting these solutions in the differential 
equation yields 

L(^> = 0 Exact solution 

-K- 

L i<t> > = R Approximate solution 
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where L is the differential operator of the problem. R is 
essentially a function whose value is equal or close to zero 
throughout the domain. Therefore this residue function is 
integrated over the domain. To pin-down the residue to small 
value at chosen number of points in the domain approximate 
weighting functions are multiplied during the minimization of 
residue in the whole solution domain. 

The above discussed objectives are acheived by setting, 

J w.R dV = 0 ...(3.2) 

D 

for i = 1 , 2 n 

where w^ are chosen weighting functions and D is the solution 
doma i n . 

Substituting equation 3.1 in 3.2 it gives 

f w. !_(<£*> dV = 0 ... <3 . 3) 

J X, 

“ D 

for i =1 , 2 n 

Equation 3.3 form the basis of all weighted residual 
methods for solving differential equations. Higher the value of n 

more accurate is the solution. The manner in which the weighting 

* 

functions are selected and the approximate solution <p is defined 
leads to many different weighted residual methods such as Ritz 
method, Galerkin method. Collocation method and the method of 
least squares. The Galerkin's method has a general applicability 
and hence been adopted here for the finite element solution 


procedure . 
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For Galerkin's method 

W. = N. ...(5.4) 

V t 

i =1 . 2 n 

where N. are shape functions of an element and n is the number of 

X 

node points in the element. 

3.3 Finite Element Procedure 

The finite element as the name indicates, implies that 
the solution domain is divided into many finite sub-domains i.e. 
elements. This discretization is advantageous particularly for 
complex domains as by considering the element as a building block 
complicated domain shapes can be represented. This is possible 
even if the elements are of any particular shape. 

The steps involved in any FEM solution are as follows : 

(1) discretization of domain into small elements 

(2) derivation of elemental properties 

(5) grouping the elements into a global assembly 
<4> incorporation of boundary conditions 

(5) solving the resulting global matrix equations to 
obtain the field variables 

(6) any processing that can further be done. 
These steps are discussed in detail in the later part 

of this chapter. The advantage of splitting the domain into many 
elements is that the field variables can be interpolated using 
standard interpolation functions within each element. For this 
purpose nodes are selected within each element and the unknown 
variable is expressed in terms of the values of field variable at 
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the selected nodes. The interpolating functions for an element 
are also called shape functions. Every node has a shape function 
whose value is unity at that particular node and zero at all the 
other nodes of that element. This is a common property of all the 
shape functions although in their form they may differ from each 
other . 

Mathematically the value of a field variable within an 
element can be expressed as. 

* T* 

4 > = 2 N iA ... <3. 5) 

i=l 

where 4 > i are nodal values of the field variable. 

Since such representations are possible for all the 
elements in the domain it can be said that the field variable is 
known in the entire domain. 

3.4 Application of FEM in the Two Phase Model 
3.4.1 Domain discretization 

Discretization of domain involves division of 
the solution domain into small elements. For plane 
two-dimensional problems, these are n-sided polygons and hence 
many types of elements are possible like triangular, 

quadrilateral etc. The domain discretization seemingly simple, 
can affect the accuracy of the solution considerably, if done 
improperly. This is more important when the domain is of complex 
shape. Isoparametric elements (so-called because the field 
variable and the spatial coordinates can be defined with the help 
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of the same set of interpolation functions) offer many advantages 
and thus used in the present investigation. Eight noded 
quadrilateral elements as shown in Fig. 0.1) have been used. The 
domain with elements is shown in Fig. <3. 2). 


3.4.2 Derivation of finite element equations 

Governing equations for the methods are given 
in chapter 2. Employing the Galerkin weighted residual method we 
get the following equations 




p<— ♦ *!“>) dA = 0 

ax 2 ay 2 J 


. . . (3.6) 



+ 


a p 

W ~ 




a 2 v 


ax 



dA = o 


. . . (3.7) 



. . . (5.8) 


where D* represents the domain of an element, ne is the number 
of elements in the domain, i varies from 1 to 8 and 1 varies from 
1 to 4. The velocity variables U and V and pressure P are 
interpolated within each element through the expressions. 


a a 

u = E N.U. , V = E N.V. and 

. V V , IV 

%=i v=± 


p = E HP, 

i=i 1 l 


. . . (3.9) 


where are quadratic interpolation functions and are linear 
interpolation functions in 2-D space. Lower order interpolation 
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Fig. 3.1(a) Eight noded element showing the local 
co-ordinate system. 




3Q 



Fig. 3.1(b) Parabolic eight noded element in the 
physical domain. 



v s. 4.0 



v = 0.0 


64 Elements 
225 Nodes 

Fig. 3.2 Finite element mesh for the solution 
domain. 
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for pressure has been done for the sake of numerical stability. 
The expressions for and are presented subsequently. The 
equations 3.6 and 3.7 are respectively usedfor the solution of 
U-velocity and V-velocity components. For the solutions of 
pressure variables, the incompressibility equation 3.8 is used 
in an indirect manner by minimizing its residue with respect to 
the pressure interpolation functions . Taken together, 
equations 5.6 to 3.8 are employed to obtain the nodal velocity 
and pressure variables. 

Once the concept of subdividing a physical domain into 
elements is accepted the problem of depicting a variation in the 
var iabl es (U , V, P> across the whole domain becomes far simpler 
since the variation can now be related within each element. Given 
such a variation within an element and the fact that elements are 
interconnected, it is a simple matter to visualise that the 
variation across the whole domain occurs in a piecewise manner. 
Interpolation functions which are generally polynomial are used 
to describe a variable within an element. The type of elements 
used to discretize the domain depends on the nature of the 
interpolation functions. For simplicity of integration and 
programming field variables as well as spatial coordinates are 
both described using interpolation functions of local 
coordinates. When same type of polynomial functions are used for 
describing both element geometry and variation of field 
variables, the element is called an isoparametric element. 
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3.4*3 Isoparametric transformation 

Any field variable can be interpolated. 

Let 

4> = oi^ + <» 2 f + n . . . ( 5 . 1 0 > 

be variation of 4> with f ,r) (local coordinates) in a four noded 
rectangular element. A natural extension of Eq.O.lO) is 


r 


r* 



> 


r. 


+3 


1 





a 

1 


+z 
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► “ 
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'z 



cx 

2 

... (3.11) 
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V * J 


<P , <p 

1 2 

, <b otnd 

~3 


are 

va lues 

of field 

variable <t> and 


<? ,7) ) , ,7) ) , ,7) ) and (£ ,7) ) are values of (? ,7?) at the 

ill 2 2 39 4 4 

four nodes of the element. 

Eq.(3.11) can be rewritten as 



* 2 

+3 



Therefore 


Cal = ter 1 i4>^ 


Mso 


r 



4> = 11 ? * tvl i 


a 

2 




<x 

4 


.. .<3.12) 


. . . <3. 13) 


. . . <3. 14) 
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Therefore <fi = C 1 ? r> ?r>] CC] _1 ...(3.15) 


or = CN N N N 3 

12 3 4 


« • • « 


where N , N , N and 

1 2 3 

N 

4 

are 

called shape 

f unct i ons 

of 

this 

particular element. 

Each 

of 

the functions 

correspond 

to 

the 

number of node given 

i n 

the l r 

subscript. They 

are equal 

to 

un i ty 


at the corresponding nodes and equal to zero at other nodes. 

For the present investigation isoparametric parabolic 
elements as shown in Fig. (3.1) has been used. 

For ease of computation, as stated earlier, the shape 
functions have been defined in terms of local coordinates. Two 
types of shape function N (quadratic) and (lower order 

interpolation) used for describing variation of variables U,V and 
P respectively are 
(a) Corner nodes: 

N. = I (1 + ^f) (?.? + 77.77 - 1) (1 + 7^7)) ...(3.17) 

Hid side nodes : 

N = j (1 - ? 2 > (1 + 7^7)), ?. = 0 

\ = j <1 + e t ?> (1 - 77 2 ), 77 . = 0 





where 'i ' is the number of nodes and , 7 ^ ) are coordinates of 
the respective nodes in the local coordinate system. 

(b) Corner nodes: 


= I (1 + (1 + 7 ^ 7 )) . . . < 5 . 1 9 ) 

where 1 is the number of corner nodes which varies from 1 to 4. 

It has been stated earlier that in the present 
investigation i soparametr i c elements have been used. 

Accordingly, 


X = P NX. 

IV 

i=l 


Y = r NY. 

I V 

1=1 


J 


. . . ( 5 . 20 ) 


where (X. ,Y. ) are spatial coordinates of the nodes. 


Comming back to Galerkin weighted residual method, using Eq.5.19 
the X-momentum equation becomes 



+ 




a*N. 

— J u. + 
ax 2 J 



= o 


. . .( 5 . 21 ) 


and Y- momentum equation becomes 
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ne r> n n n 

Vf TV V dN V Y” dN 

4 nJp z_ n k u k 4 W v j + ^jv k 4 37 J ^ 

1 J [ K=i j=i K=i j=l 


+ 


I 

1=4 


dM t 

37 



0 


... ( 3 - 22 ) 

Invoking Green.' s theorem the second order terms can again be 
reducd to the following expression : 




. . . (3.23) 


n n 

Also the terms ( £ N U £ 

K=1 j=± 


dN. n 

S)T U P and C £ 

K=i 


N V 


K K 


n dN. 

E 3>r v j3 

j=i 


are 


making the equations non-linear. So a guess value has been 
assumed for L> k and and compared with calculated values of U 
and V. This leads to an iterative procedure. Therefore X-momentum 


equation becomes : 
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and Y-momentum equation becomes 


Z f [y Z Z‘ z 


df\l 

j 
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m 



dN. dN. £N. dN. 

L L v. + I-’' —• 1 

dX dX J dY 9Y 




dN. 

N. — J V. d T 
v *n J 


J " N i 

r* 


av. 



dr 


o 


where U K and V k are guess 


values (usually taken as 


. ( 3 . 24 > 


. ( 3 . 25 ) 


prev i ous 
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iteration value or initial guess). For the continuity equation 
is used as weighting function. So 



. . . (5.26) 


The interesting point to note in equations 5.24, 5.25 and 5.26 is 
that there are twenty variables namely U and V at eight nodes and 
pressure at four corner nodes and there are twenty simultaneous 
equations per element. 

Now, the next step is to calculate the coefficients 
associated with the twenty unknowns in these equations. For 
calculating the coefficients integration of certain functions 
have to be carried out. So the concepts of numerical integration 
have been reviewed. 

3.4.4 Numerical integration 

A numerical integration procedure is adopted 
where the sampling points are used and termed as Gauss-points. In 
particular three point Gauss-Legendre quadrature has been used 
which leads to high accuracy. The unidimensional example as shown 
is Fig. <5. 5) illustrates discrete sampling points and, generally, 
the Gaussian quadrature rule leads to an equation of the form 

+ 1 m 

I = / f (?)d? « E a i f(? i > 


. . .<5.27) 
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7 6 5 


4 


1 2 3 

► 

x 

Fig. 3.3 Nodal definition of a parabolic element 
with a 3x3 Gauss point scheme. 
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where 

m 

= total number 

of integration points 


a 

V 

= ith weighting 

factor 



= coordinate of 

the ith integration point 


For two dimensional case 

+1 +1 

I = / S F(£ ,7>>d£ d7) 

-i -l 

m m 

* E E a a. F<f. , 77 . ) . . . (5.28) 

J l I V 

j=l i=± 

Clearly from equations 5.24, 5.25 and 5.26 if a local coordinate 
system is to be used then the first order variations with respect 
to the spatial coordinates must also be expressed in local terms. 
This can be acheived by chain rule. 
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in a more convenient from. 
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From equations 3.20 and 3.30 


CJ] = 


e dN. 

E 3T X v 

t=i 


8 dN. 


w~ 


a dN. 


W \ 


8 dN. 

E 75Z- Y ; 


Therefore 


detC J] 


Where 


detCJ] = 


dX dY 


dX dY 


~&fi ~Sr) ~5% 


Another transformation that is required 


integration is 


dX dY 


detCJD d£ dr? 



3.4.5 


El o merit assambl y 


SI 


In section 3.4.5 the elemental contributions 
to the left hand side coefficient matrix and right hand side 
vector were discussed. The elemental contributions are assembled 
into global matrix by adding the entries corresponding to each 
nodal variable. We get 

CA3 [X3 = CB3 . . . (3. 36> 

where [A3 = Coefficient matrix 

[X3 = Solutiuon vector 
[B3 = Right hand side vector 

Eq.(3.36) has been solved by Frontal method which has been 
described in the next section. 

There are two important features related to calculation 
of coefficients that should be noted 

<1> E anci E <Vk are actuall V obtained as guess 

in the first iteration 

(2) CA3 is calculated using 3-point gauss rule. And fj 
for each gauss point within an element is 
different . 

= K e 

The value of K (which depends on whether it lies on the matrix 
or on the second phase plate) and e (as calculated in previous 
iteration) is different at all gauss points. 

3*5 Programming 

3.5*1 Salient features of the program 

The computer program is based on the 

cental l^rary 

11 • a o; Jd 

limn iii i. 1 1 I Ill r ir«i-Ti-n*T— — t 

| dec. No. 
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algorithm developed by Taylor and HughesC13]. A very simple flow 
chart of the algorithm is depicted in Fig. 0.4). 

The main program calls the subroutines DIMENS, DINPUT, 
DRIVES, ITERAT, MATPRT and SRCAL which calls other subroutines. 

The functions of the various subroutines are as 

f ol lows : 


(a) DIMENS returns the value of the array dimensions 
which remain unchanged through out the program. 

<b> Subroutine DINPUT reads and returns all of the 
required input parameter which are listed below : 

K 

(1) slope of the second phase plate(m), L, D and 

± 

<2) ax i -symmetr i c flow indicator, the number of 
initial and boundary conditions, body forces, 
tolerance, relaxation parameters K (or K ), 

matrix 1 

dens i ty 

(5) boundary conditions. 

The subroutine DINPUT also does mesh generation by 
calling the subroutine MESANG. Also some of the given input is 
checked for correction by two routines named DIAGN1 and DIAGN2. 

<c> Then subroutine MATPRT is called. Using the values 
of m, L and D it assigns the value of variable valt at different 


gauss-points if 
phase plate it 
value one. 


the gauss-point 
K 

s assigned . 


lies within the 
Otherwise it 


domain of second 
is assigned the 


(d> After this subroutine DRIVES is called. It 
3N t dN t dM v dM 

calculates jy— and detC J] at all the 

gauss-pcjirits in all the elements. It 


further 


calls 


the 




Fig. 3.4 


Flow chart for finite element procedure 
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subroutines SHAPE4 , SHAPE 8 and DJACON. SHAPE4 and SHAPE8 

calculated M and N. at the gauss-points. DJACOB calculates the 

l dN. dN. 

values of det[J3 gg — , gj - , 

■jpg- and gg- at gauss-points. 

<e> After this the subroutine HERAT is called which 

m 

calls subroutines PRESCR, FORNTS, WRITER and TOLREL. Subroutine 

£ji 

PRESCR is called when gradient of variables such as gg ,gg are 

specified in the boundary condition. It calls another subroutine 

SURFIN. Since gradient boundary condition is not used in the 

present investigation the subroutine PRESCR has no use. 

Next subroutine that ITERAT calls is FRONTS. FRONTS 

calls subroutine MATRIX, which calculates <20 x 20) coefficient 

3N. 

matrix for a given element. It uses the values of N. , M, gy—- 
gj—i gy— . gj — . detCJ], valt and previous guess values of field 
variables (U,P,V and e which is calculated using U and V) to 
calculate coefficients. 


In FRONTS coefficient matrices for all elements are 


assembled, boundary conditions applied and matrix equation is 
finally solved using frontal method to get the value of field 
variable (U.P.V). This is discussed in details in section 

These values are then printed using subroutine WRITER. 
TOLREL compares the values of calculated field variables with the 
previous iteration values and checks whether convergence criteria 
is satisfied. If covergence criteria is not satisfied then it 
relaxes the values for the next iteration. 


<f> Finally subroutine SRCAL calculates the values of 



^ f ^ cr <y , a a , 

± 1 * 22 £2 *f • £ £1 2 2 £2 


wit 


and finds out their averages. 

The general program structure is 
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at all gauss point 

shown in Fig. (5.5). 


3*5.2 Matrix solution procedure 

Here the frontal solution method has been 
used to solve the assembled matrix equations. It is a wll known 
fact that the method adopted for solving the assembled matrix 
equation has a significant bearing on the computer storage 
requirement and execution time when the total number of unknown 
variables is very large. It becomes very time consuming and the 
core memory problem is very acute when a large number of unknowns 
are solved by the matrix inversion technique. The frontal 
solution technique can handle quite a large number of variables 
without the problem of large core memory and the computer 
execution time is also quite less. 

The frontal solution technique used in the present work 
is based on the direct Gaussian elimination procedure for solving 
the symmetric matrices where the leading diagonal is always used 
as pivot. For unsymmetric matrices encountered in a wide range of 
engineering problems, the most suitable pivot is not necessarily 
on the leading diagonal and a Frontal solution routine exists for 
off diagonal pivoting. But since this tends to be more time 
consuming, the method used in the present work uses only diagonal 
pivoting and incorporates many features of the Frontal method for 


solving symmetric matrices. 



MESANG SHAPE8I PRESCR IfRONTsI ItOLREL 


d jacobI matrix 
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Fig. 3.5 General program structure. 
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The overall Frontal solution technique covers the 
following steps : 

<i> formation of the element matrices 

(ii) assembling into a small, temporary global matrix 
of chosen size 

(iii) introduction of known variable boundary condition 

(iv> reduction of the global matrix using Gaussian 

elimination procedure 

<v> back substitution. 

The primary objective of the Frontal method is the 
elimination of variables as soon as possible after their 
introduction, via appropriate equations, into global matrix. When 
the contributions from all the elements to a particular node 
point have been assembled, the corresponding variables associated 
with that node can be eliminated. The complete matrix is 
therefore never assembled since the reduced equation can be 
eliminated from the core and stored on disc. The equations held 
in core, with the corresponding nodes and variables, are termed 
as the FRONT and the number of unknown variables in the front is 
termed as the FRONT WIDTH. 

For a non-symmetr i c global matrix, a preassigned global 
matrix core area is filled from contributing elements, the 
largest diagonal entry is the pre-ass igned core is then found and 
used as the pivot in a direct Gauss ian-e 1 imi nat i on process. As 
the maximum predetermined number of equations are eliminated, the 
corresponding reduced equations are written on to disc and more 
elements and corresponding equations are taken into core. 


The equations, nodes and variables currently in core 


are termed active, those assigned to disc 
those yet to appear in core as inactive, 
diagrammat i cal 1 y in Fig. 0.6). 


as deactivated and 
This is shown 



O Inactive Nodts 


• 

© 

Activ® Nodts I 

D®octivot*d Nod®s j 1 

V77X 

Front / j 


V//4 

Next Element of Assembly I j 

i j 


* • • 

• • • 

Assembled Element 


Fig. 3. 6 Basic idea of Frontal method 









CHAPTER 4 - 


RESULTS AND DISCUSSIONS 


4.1 Introduction 

The high temperature plastic deformation behaviour of 
two-phase alloys having the second phase with plate morphology 
has been modelled. Constitutive equation for high temperature 
plastic deformation has been taken to be that of a viscoplastic 
material, i.e. the flow stress is related to strain rate e by 
equation a = K c and pseudo viscosity of the material is given 
as /i = K e , where K and m' are material constants. The 
variable parameters of the model were , the aspect ratio, 
volume fraction and plate orientation of the second phase. 
Results of the model calculations have been generated using 
Finite Element Method. 

In this chapter the effect of various parameters during 
high temperature plastic deformation of two-phase alloys having 
plate like second phase particles on the averaged values of 
deviatoric stress components, a' and a' has been discussed. 

41 22 

Also trends of stress and strain rate fields around the second 
phase plate during hot deformation has been examined. The 
application of results of the present model to high temperature 
plastic deformation of two-phase structure are also discussed. 


4.2 Variation of Stress and Strain Components Over the Two Phase 
Domain 


The variation of a , a' , 

*q 22 


a* , 
•<* 


hydrostatic stress P, 
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and c over the two-phase domain as computed as per section 
3.5 of chapter 3 is shown in Figs .4.1 to 4.9. The figures are of 
two types i.e., 3-D plots and isolines. 

Fig.<4.1> shows the variation of a over the two-phase 
domain for K /K = 4.0, i.e. the matrix is considerably softer 

2 t 

than the second phase. L/D ratio used for generating this figure 

was kept as 1.0 with the volume fraction of the second phase kept 

as 0.24 and slope as 0.0. It can be seen that the stress 

distribution within each phase is somewhat heterogeneous but the 

extent of heterogeneity is not very high. While this minor 

heterogeneity for the softer phase is confined to the interface 

of two-phases, the heterogeneity in harder phase appears to scan 

through a much wider area fraction. Further, a is seen to be 

rising in a very steep way in the harder phase. 

Figs. 4.2a and 4.2b show isoline and 3-D plot 

respectively of the variation of the deviatoric stress component 

O’* over the two-phase domain. Similar to Fig. (4.1) average 
22 

stress is much higher for second phase plate. Also a' is seen to 

22 

be rising in a very steep way in the harder phase. 

Heterogeneity within each phase is more than that observed for 
a . Heterogeneity within the phases is in general more in harder 

•q 

phase than in softer matrix and it is maximum close to the 

two-phase boundary. The variation of a' , close to the phase 

22 

boundaries as clearly evident from Fig. <4. 2b) is somewhat 
sinuso i dal . 

Figs. 4.3a and 4.3b show isoline and 3-D plot 

respectively of the variation of the hydrostatic stress. 



0.4 


x-coo. 


i nstantaneous 
= 1 . 0 . ' 


L/D 






y-COORDINATE 


nstantaneous ft values over the domain 
1_/D =1.0, V = 0.24, m = 01 


Fig. 4. 3a Isolines for i 
(K /K = 4.0. 
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Heterogeneity in this is more than in a' and o* . Here also 

22 

close to the two-phase boundary the heterogeneity is maximum. 
Fig. 4.4 shows the variation of ex' over the two-phase domain. 

• q 

The nature of the plot is similar to that of a . 

«q 

Figs. 4.5a and 4.5b show isoline and 3-D plot 

respectively of the variation of the e . Heterogeneity in 

22 

general is much more than stress fields. Heterogeneity in the 

harder second phase is much less than that observed in the matrix 

phase. Total eight high peaks can be observed in the matrix 

phase. Four of them lie on the mid-points of the domain edges. 

Figs. 4.6a and 4.6b show isoline and 3-D plot for £ . The nature 

is same as that of e . Only a bit more heterogeneity within 

22 

second phase is observed. 

Fig. 4.7 is the isoline of a' for a domain with 
K /K = 1.5, L/D = 2.0, V, = 0.44 and plate orientation at 45°. It 

2 1 f 

is evident from this that the nature of plot for ex' is not 

•a 

affected by these model parameters. Figs. 4.8 and 4.9 are the 

isolines of e and t „ for this domain. When compared with Figs. 
22 

4.5 and 4.6 we find not much difference in nature. 

4.3 The Effect of Various Parameters of the Model on a* and ex' 

11 22 

4.3.1 Effect of change of boundary conditions 

The averaged stresses, as computed as per 
section 3.5, for a given two-phase structure (K /K , V,, L/D and 

2 1 f 

m = const.) should be invariant to the flow field or the boundary 


condition. 




















'4 


0.8 

0 . 6 ' 


0.4 







76 


Different boundary conditions (as given in appendix) 
were applied to domains with different sets of m, V , L/D values 
for K^/K^ = 4.0. The change in the value of effective pseudo 
viscosity was studied. Table 4.1 gives some representative data. 
Variation in -o' and o' with change in boundary conditions 
which give different value of average shear strain rates, is very 
small. And (u calculated using these values and and 

value show very little <~ 3%) variation. Thus it can be inferred 
that fj calculated in this way is invariant to stress and 
strain rate fields. Another interesting point to be noted from 
the table is that o' is nearly same for all the datas . p as 

• q 

calculated using o' and t cannot be said to be invariant of 


4. 3. 2 Effect of K/K 

2 1 

Effect of K /K as shown in Figs. 4.10 to 

2 1 

4.17 which is varied from two to five has been studied. The 

variation of averaged deviatoric stress components an( ^ ° f 2 .z 

with K^/K^ is almost linear. This is true for any value of m, V f 

and L/D. This fact indicates that for any morphology, the 

averaged stress-components for any two-phase can be calculated if 

the stress-component averages are known for two values of K^/K^. 

Anisotropy in general increases with K^/K^. It is substantial 

when V, is low and L/D is on higher side and m is not close to 
f 


0.5. 
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TABLE 4. 1 


Change in Effective fJ with Boundary Condition 


SI. 

No. 

K 

2 

ic - 

i 

V f 

L^D 

m 

Boundary- 

condition 

a* 

22 

£ 

22 

O'* 

•q 

• 

£ 

eq 

o> 

22 

• 

£ 

22 

O'* 

»q 

• 

£ 

eq 

1 

4 

0.44 

1.5 

0.5 

BCl 

(£ =0) 

12 

126.1 

4 

254.4 

8.0 

31.5 

31.7 

2 

4 

0.44 

1.5 

0.5 

BC2 

<£ =0.85) 

912 

122.9 

4 

255.3 

*9.2 

30.7 

27.7 

3 

4 

0.44 

1.5 

0.5 

BC3 

<£ =1.2) 

12 

122.2 

4 

255.3 

10.6 

30.5 

24.1 

4 

4 

0.44 

2.0 

0.5 

BCl 

<£ =0) 

12 

127.5 

4 

256 . 8 

8.0 

31.9 

32.0 

3 

4 

0.44 

2.0 

0.5 

BC2 

<£ =0.85) 

12 

124.2 

4 

253.1 

9.2 

31.0 

27.5 

6 

4 

0.44 

2.0 

0.5 

BC3 

U =1.2) 

12 

123.6 

4 

256.5 

10.6 

30.5 

24.2 


4. 3* 3 Effect of volume fraction 

As shown in Fig. (4. 18) for low values of 

volume fraction <<0.25>, -o' and o' vary almost linearly as a 

11 22 

function of volume fraction irrespective of morphology and 
orientation. For higher values of volume fraction -<r’ and o' 

JL % . Z Z 
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are sensitive to plate L/D and plate orientation. The variation 
of these quantities with respect to V for a given orientation 
and L/D is shown in Figs. 4.19 and 4.20. 

4.3.4 Effect of L/D 

The effect of L/D on averaged stress 
components as stated in section 2.5 is negligible for V less 
than 0.25. For higher volume fraction, as shown in Figs. 4.21 to 
4.23, averaged stresses are sensitive to L/D. The nature of 
variation is oscillating. Anisotropy is much high for high values 
of L/D except when m is equal to 0.5 <plate is oriented at 45°>. 
This effect is maximum for m = 0 (plate is parallel to domain 
edges > . 


4.3.5 Effect of plate orientation 

The effect of plate orientation is shown in 
Figs. 4.24 to 4.27. For V < 0.25 there is no significant 


variation of 

averaged deviatoric stress components 

wi th 

m. 

For 

L/D =1, -o' 

±± 

and o' variation with m is very 

22 

smal 1 

for 


varying from 

0.16 to 0.44. And for other values 

of L/D, 

-O' * 

4 4 


curve crosses the o' curve at m = 0.5. 

22 

4.4 Applications of the Results of the Model for the Analysis of 
High Temperature Deformation Processing of Two Phase Alloys 
As stated in chapter 2 the two important objectives of 
the analysis of deformation processing of two-phase alloys are : 

(a) transformation of initial microstructure to a 
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final microstructure 
(b) assessment of load requirements. 

(a) During transformation of initial microstructure to a 
final microstructure, two significant phenomenon which occur are 
shape change of second phase particle and dynamic recovery and 
recrystallization (changes taking place in matrix phase). As 
evident from Figs. 4.1 to 4.9 the deformation behaviour is 
extremely complex because of high degree of heterogeneity of 
stress and strain rate fields at microscopic level. 

The instantaneous values of o' 1 (and consequently &' ) 

22 ii 

as shown in Figs. 4.2a and 4.2b vary in somewhat sinusoidal 
manner. Since the deviatoric stress component is responsible for 
shape change this would imply that the shape of second phase 
plate like particles will evolve into a non-uniform one. 

The present model can be effectively used for a rough 
judgement of the range of heterogeneity in the stress and strain 
rate fields during the hot deformation of two-phase alloys with 
plate like second phase particle by setting the values of K^/K^ 


depending 

on 

material and trying 

out 

a 

range of V 

and 

L/D 

depending 

on 

the microstructure 

of 

the 

alloy. This 

i s 

very 

important 

i n 

connection with the 

phenomena 

of dynamic 

recovery 


and recrystallization occurring during the hot deformation of 
steel and Ti-alloys (ot+/?> . In fact temperature of deformation 
processing on which the present model parameters like K^/K^, 
and L/D depend to a large extent can give significant 
informations using the present model. From the results of the 
present model, it can be seen that on the average strain rate in 

1 
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the matrix phase is much higher than that in second phase plate 
and the range of heterogeneity of strain rate field in the second 
plates is much lower compared to that in the matrix phase. This 
is true for both K^/K^ as high as 4 (Figs. 4.5 and 4.6) and as 
low as 1.5 (Figs. 4.8 and 4.9) which is close to the value of 
in alloys like Ti-6A1-4V. This can be probable cause of the 
fact that during the thermomechanical processing of coarse 
two-phase alloys, one of the phase softens by both dynamic 
recovery and recrystallization whereas the other phase undergoes 
only recovery[8]. 

(b) It has been discussed in chapter 2 that the knowledge 
of material properties used in constitutive equation of the two 
phase material is an essential requirement to assess the load 
requirements. These material properties as stated in chapter 1 
are sensitive to the mi crostructural parameters like second phase 
morphology which changes with the deformation, making the 
assessment of load requirement very difficult. So it is very 
important to correlate the material properties used in the 
constitutive equation (or in governing equation) to the 
mi crostructural parameters like plate size and shapes. 

As stated in chapter 2 microstructure of the two-phase 
alloy as shown in Fig. (2.1a) is similar to the composite shown in 
Fig. (2.1b). Since the constitutive equation for a given two-phase 

structure (i.e. for fixed K /K ratio, L/D ratio, volume fraction 

2 1 

and plate orientation) must remain same for different flow 
fields, it must be expected that a given homogenization scheme 
must yield averaged stress values, as per computed by section 



00 

5.5, which are invariant to the flow field and boundary 

conditions, (which is the material property for viscoplastic 

flow) calculated using o' and o' as discussed in section 4.3.1 

11 22 

has been found to be invariant of flow and the u _ thus 

eff 

calculated using the present model for different zones, can be 
used to calculate the ju for the microstructure in Fig. (2. lb) 
by using the FEM program described in section 3.5 with minor 
modifications. 



CHAPTER 5 


CONCLUSIONS 

High temperature plastic deformation behaviour of 
two-phase alloys having the second phase with plate morphology 
has been modelled using Finite Element Method (FEM). Constitutive 
equation for high temperature plastic deformation has been taken 

to be that of a viscoplastic material, i.e. the flow stress a is 

✓ 

_ m 

related to strain-rate e by the equation a = K e and pseudo 

m— l 

viscosity of the material p is given as = K “e , where K and m 
are material constants . The variable parameters of the model were 
K /« , the aspect ratio, volume fraction and plate orientation of 
the second phase. Boundary condition has been also varied. The 
following conclusions can be drawn from the results of the 
present study r 

<1> Plots of various stress and strain-rate components 
over the domain, show a high level of heterogeneity which exists 
in the vicinity of second phase plate like particles. 

(2) Heterogeneity is more for the strain-rate 
components. On the average strain-rate is much higher in the 
softer matrix phase. This finding can perhaps explain the fact 
that during hot deformation of two-phase alloys dynamic recovery 
is observed in both the phases whereas dynamic recrystal 1 i zat i on 
is restricted to only one phase. 

(3) The nature of variation in the deviatoric stress 
component shows a somewhat sinusoidal nature at the two-phase 
boundary. This indicates that the shape change in the second 
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phase plate like particles will be non-uniform. 

<4> Pseudo-viscosity, /j , in the viscoplastic 

deformation model, calculated using and t , can be 

22 11 11 22 

treated as the effective pseudo-v i soc i ty , /j , of the two-phase 

structure. This fact can be utilized for calculating p ^ for a 

real microstructures in two-phase alloys. The assessment of load 

requirement for carrying out the plastic deformation of two-phase 

alloys can thus be made from the mi crost ructura 1 information. 

(5) Results showing the variations in a' and cr' with 

11 22 

model parameter indicate the following, 

<i> if fjt for a given morphology is known for 

two values of K /K the u _ for any 

2 1 eff 

two-phase structure with same morphology and 

different K /K can be conveniently 

2 1 

i nt erpolated 

<ii> for V < 0.Z5, the effect of morphology of 
second phase on fj is negligible. 
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Appendix A 



a 0.5 vf = . 44 l/d=l 5 

0 34 8 8 1 l 9o.o 

0 86 0 0 

1 -01 100 0 0 

°-* *4062 .2708 4.0 

~~ 3 4 


209 

209 

210 
210 
211 
211 
212 
212 1 
213 3 

213 

214 

214 

215 

215 

216 
216 1 
217 3 

217 1 

218 3 
218 1 
219 

219 

220 
220 
221 
221 
222 
222 


223 3 


4.0 


223 1 4 

224 3 4 

224 1 4 

225 3 4 
225 1 4 

1 3 0 

2 3 0 

3 3 0 

4 3 0 

5 3 0 

6 3 0 

7 3 0 

8 3 0 

9 3 0 

10 3 0 

11 3 0 

12 3 0 

13 3 0 

14 3 0 

15 3 0 

16 3 0 

17 3 0 
1 1 2 

18 1 2 
27 1 2 
44 1 2 
53 1 2 
70 1 2 
79 1 2 
96 1 2 



105 1 2 
122 1 2 
151 1 2 
148 1 2 
157 1 2 

174 1 2 
185 1 2 
200 1 2 
209 1 2 
17 1 -2 
26 1 -2 
45 1 -2 
52 1 -2 
69 1 -2 
78 1 -2 
95 1 -2 
104 1 -2 
121 1 -2 
150 1 -2 
147 1 -2 
156 1 -2 

175 1 -2 
182 1 -2 
199 1 -2 
208 1 -2 
225 1 -2 
1 2 0.0 



n ?A 5 n Vf= * 44 4.0 

0 39 8 8 1 1 90.0 

0 102 0 0 

1 .3 .01 1 100 0 0 
0.3 .4062 .2708 4.0 
209 3 4 


209 1 

210 3 
210 1 
211 ? 
211 1 
212 3 
212 1 
213 3 

213 1 

214 3 
214 1 
213 3 
213 1 
216 3 
216 1 
217 

217 

218 
218 
219 

219 

220 
220 
221 
221 1 
222 3 
222 1 
223 3 

223 1 

224 3 
224 1 
223 3 
223 1 
1 3 

3 
3 
3 
3 
3 
3 
3 
3 


2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 

4 

2 


0 

0 

0 

0 

0 

0 

0 

0 

0 


10 3 

11 3 

12 

13 

14 
13 
16 
17 
1 1 


18 

27 

44 

33 

70 

79 

96 


103 1 2 
122 1 2 



131 1 2 
148 1 2 
137 1 2 
174 1 2 
183 1 2 
200 1 2 
209 1 2 
17 1 -2 
17 3 2 

26 1 -2 
26 3 2 
43 1 -2 
43 3 2 
32 1 -2 
32 3 2 
69 1 ~2 
69 3 2 
78 1 -2 
78 3 2 
93 1 -2 
93 3 2 
104 1 -2 
104 3 2 
121 1 -2 
121 3 2 
130 1 -2 
1 30 3 2 
147 1 -2 
147 3 2 
136 1 -2 
1 36 3 2 
173 1 -2 
173 3 2 
182 1 -2 
1 82 3 2 
199 1 -2 
199 3 2 
208 1 -2 
208 3 2 
223 1 -2 
1 2 0.0 



a 0.5 vf= .44 1 /d=l. 5 4.0 
0 54 8 8 1 1 90.0 

0 B6 0 0 

1 .5 .01 100 0 0 

0 .5 .4062 .2708 4.0 
209 3 4 

209 1 5 

210 3 4 

210 1 5 

211 5 4 

211 1 5 

212 3 4 
212 1 5 
215 3 4 
215 1 5 
214 3 4 

214 1 5 

215 3 4 

215 1 5 

216 3 4 

216 1 5 

217 3 4 

217 1 5 

218 3 4 

218 15 

219 3 4 

219 15 

220 3 4 

220 1 5 

221 3 4 
221 1 5 
22 2 3 4 
22 2 1 5 
22 3 3 4 
22 3 1 5 
22 4 3 4 
22 4 1 5 
22 5 3 4 
22 5 1 5 
1 3 0 

2 3 0 

3 3 0 

4 3 0 

5 3 0 

6 3 0 

7 3 0 

8 3 0 

9 3 0 

10 3 0 

11 3 0 

12 3 0 

13 3 0 

14 3 0 

15 3 0 

16 3 0 

17 3 0 
1 1 2 

18 1 2 
27 1 2 
44 1 2 
53 1 2 
70 1 2 
79 1 2 
96 1 2 
105 1 2 



122 1 2 
131 1 2 
148 1 2 
137 1 2 
174 1 2 
183 1 2 
200 1 2 
209 1 2 
17 1 -2 
26 1 -2 
43 1 -2 
32 1 -2 
69 1 -2 
78 1 -2 
93 1 -2 
104 1 -2 
121 1 -2 
130 1 -2 
147 1 -2 
156 1 -2 
173 1 -2 
182 1 -2 
199 1 -2 
208 1 -2 
225 1 -2 
1 2 0.0 





implicit real *8 (a-h , o-z ) 

DIMENSION BOUDV< 624 ) , COORD ( 264 , 2 ) , EQRHS(624 > , GFLUM< 105,105), 
GRADB! 264,4) ,LB0UD(624) .LHEDVU05) .LNODS (73, 8) , 
NADFM< 264 > .NGRAD! 73) , NODFM! 264) ,PNORM(105) , 

POSGP< 3 ) , VARB1 ( 624 ) , VARB2 ( 624 ) , WE IGP( 3 > 

.STRll (73) , STR22(73 > , STR1 2 (73 ) , VALT (73 .9) 

OPEN (UN I T= 7 , F ILE= ' out 8 ' ) 

OPEN (UN I T=40 , F ILE = * i 80 . inp* ) 

OPEN (UN I T— 1 0 , F ILE= 'vall3.dat' , f orm= * unformatted * > 

OPEN(UNIT=2? , F I LE= * val 1 6 . dat ' , f orm= * un format t ed * > 

open (uni t = 8, f i le= ’sig’ ) 

open(uni t = 9 , f i le=’sref ’ ) 

open(umt=ll , f i le = ’tao’ ) 

open(un 1 1= 12 , f i le= ’ taol 1 * > 

open (uni t = 15 , f i Ie = * tao22 ’ ) 

open (uni t=14,file= , pres’> 

open (un it=15,file='srll*> 

open (uni t=16 ,fiie= , sr22*> 

open(unit=17,f i le='srl2' > 

open (uni t= 18 , f i 1 e= ' taol 2 * ) 

SET DYNAMIC DIMENSION VALUES 

CALL D I MEN S ( MELEM , MFRON , MPO IN, MTOTV) 

READ IN ALL PROBLEM DATA 


CALL D I NPUT ( BOUDV , COORD , DENSY . GRADB , I AXSY , LBOUD , LNODS , 
MELEM , MPO I N , MTOTV , NADFM , NBCON , NDOFM , NELEM , 
NEVAB , NGAUS , NGRAD , N I TER , NNODL , NNODP . NODFM , 
NPO I N , NTOTV , RELAX , TOLER , VARB 1 , VARB2 , V I SCY , 
XFORC , YFORC , SLOPE 1 . D . E , FAC . XMAX , YMAX ) 


CALCULATE THE SHAPE FUNCTION AND DERIVATIVE VALUES 

CALL DRIVES! COORD , LNODS , MELEM , MPO I N , NELEM , NGAUS , NNODL . 
NNODP , POSGP , WE I GP ) 

CALL MATPRT ( COORD , VALT . LNODS , MELEM , MPO I N , NELEM , NGAUS , 
POSGP , XMAX , YMAX , SLOPE 1 , D , E , FAC , NNODP > 

SET UP EQUATIONS AND ITERATE UNTIL SOLUTIONS CONVERGE 

CALL ITERAT(BOUDV, COORD, DENSY,EQRHS,GFLUM, GRADB, IAXSY, 
LBOUD , LHEDV , LNODS . MELEM , MFRON , MPO I N , MTOTV , 
NADFM . NBCON , NDOFM , NELEM , NEVAB . NGAUS . NGRAD , 
NITER, NNODL , NNODP , NODFM , NPO I N . NTOTV , PNORM , 
POSGP , RELAX , TOLER , VARB 1 , VARB2 , V I SCY , WE I GP . 
XFORC, YFORC, VALT) 
write!*-,*) 'no* 

CALL SRCAL < VARB 1 , LNODS , NADFM , NODFM . MTOTV , MPO I N . MELEM , 
NPO I N , NELEM , NNODP , NNODL , STR 1 1 , STR 2 2 , STR 1 2 , V I SCY , VALT 
, COORD, posgp> 



oonoo non 


wr i te (* , *> 'yes' 

CLOSE (UNI T=40> 

CLOSE (UNI T=1 0 ) 

CLOSE(UNIT=25 > 

CLOSE (UNI T= 7 > 
ciose(unit=8> 
ciose(um t = 9) 
close(uni t=l 1 ) 
c iose(un 1 1 = 12 > 
close(unit=13> 
close(unit=14> 
close(uni t = l 5> 
close(unit=16) 
close(unit=17> 
c lose (un i t = l 8 ) 

STOP 

END 

SUBROUTINE 0 1 MENS ( MELEM , MFRON , MPO I N , MTOTV > 
MELEM=73 
MFRON=l 05 
MPOIN=264 
MTOTV=624 
RETURN 
■ END 

SUBROUT I NE SHAPE4 ( DER I V , SHAPE , XEQ I V , YEQ I V ) 
DIMENSION DER I V ( 2 , 8 ) , SHAPE ( 8 > 

LINEAR ELEMENT ANTI -CLOCKWISE NODE NUMBERING 

X=XEQI V 
Y=YEQI V 
XY=X*Y 

SHAPE ( 1 ) = ( 1 -X- Y+XY >*0.25 
SHAPE ( 2 > = (1 +X- Y-XY > *0 . 2 5 
SHAPE ( 5 > = ( 1 +X+ Y+XY >*0.25 
SHAPE ( 4 > = ( 1 -X+ Y-XY > *0 . 2 5 
DERIV(1 ,1>=(-1+Y>*0.25 
DERIVd ,2>=(1-Y>*0.25 
DERIVd ,?) = (1+Y>*0. 25 
DERIVd ,4) = (-l-Y>*0.25 
DERIV(2,1>=(-1+X>*0.25 
DERIV(2,2)=(-1-X>*0.25 
DER IV(2,5>=(1+X>*0.25 
DERIVd , 4 > = (1 -X>*0 . 25 
RETURN 
END 


SUBROUT I NE SHAPES ( DER I V , SHAPE , XEQ I V , YEQ I V > 
DIMENSION DER I V ( 2 , 8 ) , SHAPE ( 8 > 

PARABOLIC ELEMENT ANTI -CLOCKWISE NODE NUMBERING 


2 



n aa ana 


c 


X=XEQI V 

Y=YEQI V 

XY-X*Y 

XX=X*X 

YY=Y*Y 

XXY=XX*Y 

XYY=X*YY 

X2=X*2 

Y2=Y*2 

XY2=XY*2 

SHAPE ( 1 > = < -1 +XY+XX+YY-XXY-XYY >*0.25 

SHAPE ( 2 > = ( 1 -Y-XX+XXY >*0 . 5 

SHAPE ( 5 > = ( -1 -XY+XX+YY-XXY+XYY > *0 .25 

SHAPE < 4 > * < 1 +X-YY-XYY >*0.5- 

SHAPE ( 5 > = C - 1 +XY+XX+YY+XXY+XYY >*0.25 

SHAPE < 6 > = < 1 + Y-XX -XXY >*0.5 

SHAPE C 7 > = ( - 1 -XY+XX+YY+XXY-XYY >*0.25 

SHAPE ( 8 > = < 1 -X-YY+XYY >*0.5 

DER I V( 1 . 1 >=< Y+X2-XY2-YY>*0 . 25 

DERI V< 1 , 2 > = < -X+XY > 

DERIV<1 , 3 > = < -Y+X2-XY2+YY>*0 . 25 

DER IV(1 , 4 > = ( 1 - YY >*0.5 

DERI V( 1 , 5 > = C Y+X2+XY2+YY>*0 . 25 

DERI V( 1 , 6 > =-X-XY 

DERI V( 1 , 7 > = ( -Y+X2+XY2-YY)*0 . 25 

DER IV(1 ,8> = (- l +Y Y >*0.5 

DER I V < 2 . 1 > = ( X+Y2-XX-XY2 >*0.25 

DER IV(2,2> = C-1 +XX >*0.5 

DERI V( 2 , 3 > = < -X+Y2-XX+XY2 >*0 . 25 

DER I V < 2 . 4 > = - Y-XY 

DER I V < 2 , 5 > = ( X+Y2 +XX+XY2 >*0.25 

DERI V< 2 , 6 > = < 1 -XX>*0 . 5 

DERI V< 2 , 7> = < -X+Y2+XX-XY2 >*0 . 25 

DER I V< 2 , 8 > = -Y+XY 

RETURN 

END 


SUBROUTINE D JACOB (COORD, DER I V, DET JB, DJ AC I ,D JACK, IELEM.LNODS, 
1 MELEM , MPO I N , NNODP > 

DIMENSION COORDCMPOIN , 2 > , DERI V( 2 , 8 > , D3ACI < 2 , 2 > , DJACK( 2 , 2 > , 

I LNODS < MELEM , 8 > 

SET UP TEMPORARY MATRIX TO ALLOW JACOBIAN TO BE FORMED 

DO 20 IDIME=I , 2 
DO 20 JDIME=1 , 2 
TEMPY=0 

DO 10 INODP=l .NNODP 
KPOIN=IABS<LNODS( IELEM, INODP>> 

TEMP Y=TEMP Y+DER I V ( I D I ME , I NODP > *COORD ( KPO I N . JD I ME > 

10 CONTINUE 

D J ACK < I D I ME . JD I ME > = TEMP Y 



20 CONTINUE 

DET JB=DJACK< 1 , 1 ) *0 JACK < 2 , 2 > -D JACK < 2 , 1 >*DJACK<1 .2) 

C 

C CHECK FOR ZERO OR NEGATIVE DETERMINANT 

C 

IF(DETJB>30, 30,40 
30 CONTINUE 

WRITE<7 , 2000 ) IELEM 
STOP 

40 CONTINUE 

C 

C INVERT TEMPORARY MATRIX TO FORM JACOB I AN- INVERSE 

C 

DJACI < 1 , 1 > =0 JACK (2,2) /DET JB 
DJACI (2,2) =0 JACK < 1 , 1 >/DETJB 
DJACI <1 ,2>=-OJACK(l , 2>/DETJB 
DJACI ( 2 , 1 > =-DJACK( 2 , 1 > /DET JB 

2000 FORMAT < / /37H NON POSITIVE DETERMINANT FOR ELEMENT, I4> 

RETURN 

END 

C 

C 

C 

SUBROUT I NE DR I VES < COORD . LNODS , MELEM , MPO I N , NELEM , NGAUS , NNODL , 
1 NNODP , POSGP , WE I GP > 

C 

C EXTERNAL SUBROUTINES: 

C SHAPE8.D JACOB, SHAPE4 

C 

DIMENSION AREAWC 9 ) ,CARLG(2 , 36 > ,CARPG< 2 , 72 ) , CARTP< 2 , 8 > , 

1 DERIV(2,8> .DJACI (2,2) ,DJACK( 2,2) ,POSGP<3> , COORD (MPO IN, 2 ) , 

2 SHALG( 36 > ,SHAPE<8) , SHAPG (72>,WEIGP<3>, LNODS < MELEM , 8 ) 

C 

C REWIND TAPE BEFORE WRITING ON SHAPE FUNCTIONS 

C 

REWIND 10 

C 

C SET UP POSITIONS AND WEIGHTS FOR 3 POINT GAUSS RULE 

C 

P0SGP<1 )=0 . 7745966692 

POSGP < 2 ) = 0 

POSGP < 3 > = -POSGP ( 1 > 

WEIGP<1)=0. 5555555556 
WEIGP( 2 >=0.8888888889 
WE I GP ( 3 > =WE I GP ( 1 > 

C 

C CALCULATE SHAPE FUNCTIONS AND LOCAL DERIVATIVES FOR ELEMENTS 

C 

DO 60 IELEM=1 .NELEM 
LGAUS=0 

DO 50 IGAUS=1 .NGAUS 
DO 50 JGAUS=1 .NGAUS 
LGAUS=LGAUS+1 
XEQIV=PQSGP( IGAUS) 

YEQ I V=POSGP ( JGAUS > 


4 



c 

C USE GAUSS POSITIONS TO CALCULATE LOCAL VALUES 

C 

CALL SHAPES < DER I V , SHAPE , XEQI V , YEQIV) 

C 

C SET UP JACOBIAN MATRIX AND INVERSE 

C 

IELEMS= IELEM 

CALL DJ ACOB< COORD .DERI V , DET JB , D JACI , DJACK , IELEMS , LNODS , 
1 MELEM , MPO I N , NNODP ) 

C 

C CALCULATE GLOBAL DERIVATIVES AND AREA*GAUSS WEIGHTS 

C 

DO 10 IDIME=1 , 2 
DO 10 I NODP = 1 , NNODP 
CARTP< IDIME , INODP) =0 
DO 10 JOIME=l , 2 

CARTP( IDIME . I NODP >=CARTP (IDIME, INODP >+ 

1 DJACI (IDIME , JDIME>*DERIV( JDIME , INODP) 

10 CONTINUE 

ARE AW ( LGAUS > =DETJB*WE I GP ( I GAUS ) *WE I GP ( JGAUS > 

C 

C PUT SHAPE FUNCTIONS AND DERIVATIVES IN ELEMENT MATRIX 

C 

DO 20 I NODP = 1 , NNODP 

KGAPA = ( LGAUS - 1 > *NNODP+ 1 NODP 

SHAPG ( KGAPA ) =SHAPE ( I NODP ) 

DO 20 IDIME=1 . 2 

CARPG (IDIME. KGAPA > =CARTP ( I D I ME . I NODP ) 

20 CONTINUE 

C 

C USE GAUSS POSITIONS TO CALCULATE LOCAL VALUES 

C 

CALL SHAPE4 ( DER I V , SHAPE , XEQI V. YEQIV) 

C 

C CALCULATE GLOBAL DERIVATIVES FOR LINEAR FUNCTIONS 

C 

DO 50 IDIME=1 , 2 
DO 50 INODL= 1 , NNODL 
CARTP( IDIME, INODL>=0 
DO 50 JDIME=1 , 2 

CARTP( IDIME. INODL>=CARTP( IDIME, INODL>+ 

1 DJACI ( IDIME, JDIME>*DERIV( JDIME, INODL) 

50 CONTINUE 

C 

C PUT SHAPE FUNCTIONS ANO DERIVATIVES IN ELEMENT MATRIX 

C 

DO 40 I NODL =1 , NNODL 

KGAL I = ( LGAUS - 1 >*NNODL+ INODL 

SHALG ( KGAL I > = SH APE ( I NODL ) 

DO 40 IDIME=1 , 2 

CARLG (IDIME, KGAL I ) = CARTP ( I D I ME . I NODL ) 

40 CONTINUE 

50 CONTINUE 

C 
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C NOW WRITE ALL THE ELEMENT SHAPE FUNCTIONS ON TAPE 

C 

NGAPA =LGAUS *NNODP 
NGAL I =LGAUS*NNODL 

WR I TE { 1 0 ) I ELEM , LGAUS , NGAPA , NGAL I . 

1 ( (CARPG( IDIME. IGAPA) , IDIME= 1 . 2 > , SHAPG( IGAPA) , IGAPA=1 .NGAPA) , 

2 ( ( CARLG ( I D I ME , I GAL I > , IDIME = 1 ,2) , SHALG( I GAL I ) . I GAL 1*1, NGAL I > . 

3 ( AREAW < I GAUS > , IGAUS = 1 .LGAUS) 

60 CONTINUE 

RETURN 
END 
C 
C 
C 
C 
C 

SUBROUT I NE SRC AL ( VARB 1 , LNODS , NADFM . NODFM , MTOTV , MPO I N . MELEM , 
1 NPO I N , NELEM . NNODP . NNODL . STR 1 1 , STR2 2 , STR 1 2 , V I SCY , V ALT , COORD 
1 .posgp) 

DIMENSION STR1KMELEM) . STR22 (MELEM) , STR1 2 <MELEM> 

1 , CARTP( 2 , 8 > , NODFM(MPOIN) ,NADFM<MPOIN> . LNODS (MELEM , 8 ) 

1 , VARB1 (MTOTV) , SHAPG< 72 ) , CARPGC2 , 72 ) , UVEL(8 > , VVEL( 8 > , 

1 SHAPL ( 4 ) . SHAPP ( 8 ) , SHALG( 36 ) , CARTL (2,4), CARLG (2.36) 

1 , AREAW ( 9 ) , VALT (MELEM, 9 ) , COORD ( MPO I N , 2),srll(melem,9) , 

1 sr22(mel em, 9>,srl2(melem,9),tll (me lem, 9),teq(melem,9>, 

1 t22(melem,9),tl2(melem,9),press(melem,9),sref(melem,9>, 

1 sigll(melem,9),sig22(melem,9),sigl2(melem,9).sigeq(melem,9> 
1 ,deriv(2,8) , shapx(8> ,posgp(3> , gcoordl (me 1 em, 9 ) . 

1 gcoord2 (me 1 em, 9) ,areawe(melem,9> 

REWIND ( 10) 

DO 110 IELEM=1 .NELEM 
STR1 1 ( IELEM) =0 
STR22(IELEM)=0 
STR 12(1 ELEM ) = 0 

READ (10) JELEM . LGAUS . NGAPA , NGAL I , 

1 ( (CARPG( IDIME , IGAPA) , IDIME=1 , 2 ) , SHAPG( IGAPA) , IGAPA=1 , NGAPA) , 

2 ( ( CARLG ( IDIME , I GAL I ) , IDIME=1 , 2 ) . SHALG( I GAL I ) , IGALI = 1 .NGAL I ) . 

3 ( AREAW ( I GAUS ) , IGAUS=1 .LGAUS) 
write(*,*> * lgaus= ' , lgaus 
DO 100 IGAUS=1 .LGAUS 
DAREA=AREAW( IGAUS) 
areawe ( i e 1 em , i gaus > =darea 
DO 30 INODP=l , NNODP 

SHAPP ( I NODP > = SHAPG ( NNODP* ( I GAUS - 1 ) + 1 NODP ) 

DO 30 IDIME=1 , 2 

CARTP ( I D IME , I NODP ) =CARPG (IDIME. NNODP* ( I GAUS - 1 ) + I NODP ) 

30 CONTINUE 

DO 40 INODL=l , NNODL 

SHAPL ( I NODL ) = SHALG ( NNODL* ( I GAUS - 1 ) + I NODL ) 

DO 40 IDIME=1 , 2 

CARTL ( I D IME , I NODL ) =CARLG (IDIME. NNODL* ( I GAUS - 1 ) + I NODL ) 

40 CONTINUE 

UVELY=0 
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VVELY=0 

srll( ielem, i gaus ) = 0 
sr22 < i e 1 em, i gaus ) =0 
srl2< ielem, igaus>=0 
sref ( ielem, igaus)=0 
C 

do 41 i nodp=l , nnodp 

KPO I N= I ABS ( LNODS ( I ELEM , I NODP > ) 

I TOTU=NADFM< KPO I N ) 

I TOTV= I TOTU+NODFM < KPO I N > - 1 
CARX=CARTP(1 , I NODP) 

CARY=CARTP (2,1 NODP > 

srll< ielem, i gaus ) =varbl < i totu >*carx + srl 1 ( i e 1 em, igaus > 
sr22< ielem, i gaus ) =varbl ( i totv )*cary+sr22 ( ielem, igaus > 
srl2< ielem, i gaus )=varbl < i totu)*cary*. 5+varbl ( i totv)*carx* . 5 
1 +srl2( ielem, igaus) 

41 continue 

srllav=srl lav+srl 1 ( i elem, i gaus )*darea 
sr22av=sr22av+sr22< ielem, igaus)*darea 
srl2av=srl2av+srl2( ielem, igaus)*darea 
sref< ielem, i gaus ) = sqrt (2*(srll(ielem, igaus )**2 ) 

1 +(srl2< ielem, igaus)**2) 

1 +2*( sr22 ( i e 1 em, i gaus)**2 ) ) 

sref av=s ref avFs ref ( ielem, igaus)*darea 
vsfac=sref( ielem, l gaus )**- . 95 

tl 1 ( i e lem, i gaus ) =vsfac*val t ( i e 1 em, i gaus )*v i scy* 

1 srl 1 ( i e 1 em, i gaus ) 

t22 ( ielem, igaus ) =vsf ac*val t ( i e lem, l gaus )*v i scy* 

1 sr22( ielem, igaus) 

tl2< ielem, igaus ) =vsfac*val t ( ielem, i gaus )*v i scy* 

1 srl2( ielem, igaus) 

teq< ielem, igaus)=valt( ielem, igaus)*viscy*vsfac* 

1 sref < ielem, igaus) 

teqav=teqav+teq< ielem, igaus)*darea 
C EVALUATE VELOCITIES & STRAIN-RATES AT GAUSS POINTS 


c EVALUATE HYDROSTATIC STRESS AT GAUSS-POINTS 

press < ielem, i gaus ) = 0 . 0 
DO 49 I NODL = 1 , NNODL 
JNODL =2*1 NODL - 1 

KPO I N= I ABS < LNODS < I ELEM , JNODL ) ) 

I TOTU = NADFM < KPO I N ) 

ITOTP=ITOTU+l 

press ( ielem, igaus ) =press( ielem, i gaus ) 

1 +varbl ( i totp)*shapl ( inodl ) 

49 CONTINUE 

pressav=pressav+press( ielem, igaus)*darea 
C 

DO 50 INODP=l , NNODP 

KPO I N= I ABS < LNODS ( I ELEM , I NODP > ) 

I TOTU = NADFM ( KPO I N ) 

I TOT V= I TOTU+NODFM ( KPO I N ) - 1 
SHAPE=SHAPP ( I NODP) 

CARX=CARTP( 1 . I NODP) 
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CAR Y=CARTP (2,1 NODP > 

UVELY=UVELY+VARB1 ( ITOTU)*SHAPE 
WELY=WELY+VARB1 ( ITOTV)*SHAPE 

STR1 1 ( IELEM) =STR1 1 C IELEM) 4-VARB1 ( ITOTU)*CARX*DAREA*VISCY 
1 *VALT ( IELEM, IGAUS)*vsfac 

STR22( IELEM) =STR22( I ELEM) 4-VARB1 < ITOTV)*CARY*DAREA*VISCY 
1 *VALT < I ELEM , I GAUS ) *vs f ac 

STR1 2 ( I ELEM) =STR1 2 ( I ELEM) 4-0 . 5*( VARB1 ( I TOTU ) *CARY 
1 4-VARB1 ( ITOTV)*CARX)*DAREA*VISCY*VALT( IELEM, I GAUS )*vsfac 
50 continue 

100 continue 

110 continue 

do 112 ielem=l .nelem 
do 111 kgaus=l,9 

s igl 1 < i el em, kgaus ) = -press ( i e 1 em, kgaus )4-t 1 1 ( i e 1 em, kgaus ) 
s ig22( i e 1 em, kgaus ) = -press < ielem,kgaus)4-t22( ielem, kgaus) 
sigl2< ielem, kgaus )=t 12 ( ielem, kgaus) 

s i geq( i el em, kgaus )=sqrt((sigll(iel em, kgaus) -s ig22 < i e lem 
1 , kgaus ) )**2 ,04-6.0*(sigl2( ielem, kgaus >**2 .0) 

1 +s igl 1 < i el em, igaus )**2 . 0 ) 

s i geqav = s igeqav4-s l geq( ielem, i gaus )*areawe( ielem, i gaus ) 

111 continue 

112 continue 

do 115 ielem=l , nelem 
lgaus=0 

do 114 igaus=l,5 
do 115 j gaus =1,5 
xeqiv=posgp< igaus) 
yeqiv=posgp( jgaus) 
lgaus = lgaus4-l 

call shape8(der i v , shapx , xeqi v , yeq i v ) 
gcoordl (ielem, Igaus ) = 0 . 0 
gcoord2( ielem, lgaus)=0 . 0 
do 1150 inodp=l .nnodp 
kpoin=iabs< lnods( ielem, inodp) ) 
write (70,*) kpoin 

gcoordl (ielem, Igaus ) =gcoordl (ielem, Igaus ) 4-coord ( kpo in , 1 > 

1 *shapx( i nodp) 

gcoord2( ielem, Igaus) =gcoord2 (ielem, Igaus) 4-coord (kpo i n , 2 ) 

1 *shapx( i nodp) 

write(70,*) shapx ( i nodp) , coord (kpo l n , 1 ) , coord (kpo in , 2 ) 

1150 continue 

115 continue 

114 continue 

115 continue 
avl 1 =0 
av22=0 
avl 2=0 

DO 120 IELEM=1 .NELEM 

WR I TE ( 7 , * ) I ELEM , STR 1 1 ( I ELEM ) , STR22 ( I ELEM) , STR 1 2 ( I ELEM ) 
avl l=avl 14-strl 1 ( ie 1 em) 
av22=av224-str22( Ielem) 
avl2=avl24-strl2( ielem) 

120 CONTINUE 

WR I TE ( 7 , * > srl lav , sr22av , srl 2av . sref av , teqav ,avl 1 , av22.avi 2 
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WRITE(*,*> srllav, sr22av , srl 2av , sref av , teqav.avl 1 ,av22 ,avl2 
do 125 i e 1 em=l , ne lem 
do 124 igaus = l , 9 

c wr i te<8 , *) gcoordl ( i e 1 em, i gaus ) , gcoord2 < i e 1 em, i gaus ) 

c 1 ,sigeq< ielem, igaus) 

c write<9,*) gcoordl ( i e lem, i gaus) , gcoord2 ( i e 1 em, i gaus > 

c 1 , sref ( i e 1 em, i gaus ) 

c writedl,*) gcoordl ( i e lem, i gaus ), gcoord2 ( i e 1 em, igaus ) 

c 1 , teq< ielem, igaus) 

c write(12,*> gcoordl < i e 1 em, i gaus ), gcoord2 ( i e 1 em, i gaus ) 

c 1 , 1 1 1 ( i e 1 em, i gaus > 

c write<13,*> gcoordl ( ielem, igaus) ,gcoord2( ielem, igaus) 

c 1 , t22< ielem, igaus) 

c write<14,*) gcoordl ( i e 1 em, i gaus ), gcoord2 ( i e lem, igaus ) 

c 1 ,press( ielem, igaus) 

c write<15,*) gcoordl < l e 1 em, i gaus ), gcoord2( i e 1 em. igaus ) 

c 1 , srl 1 < i e 1 em, i gaus ) 

c write(16,*> gcoordl ( i e 1 em, igaus ), gcoord2< i e lem, i gaus ) 

c 1 , sr22< i el em, igaus ) 

c write<17,*) gcoordl ( i e lem, i gaus ), gcoord2 ( i e lem, igaus ) 

c 1 , srl 2 < i el em, i gaus ) 

c wr i te ( 1 8 , *) gcoordl ( ielem, igaus) ,gcoord2( ielem, igaus) 

c 1 , tl 2 ( i e lem, i gaus ) 

124 continue 

125 continue 

c write (8,*) sigeqav 

write(9,*) srefav 
writedl ,*) teqav 
write(15,*) srllav 
wri te<16 , *) sr22av 
RETURN 
END 


SUBROUTINE SURF IN(COORD , COSLX , COSLY , I ELEM, IGAUS . IS IDE ,LNODS , 

1 MELEM , MPO I N , NNODP , POSGP , SHAPE , SLETH , WE I GP > 

EXTERNAL SUBROUTINES : 

SHAPE8.D JACOB 

DIMENSION DERI V< 2, 8) ,D JACK 2, 2) ,0 JACK< 2 , 2) , POSGP< 5 > , WEIGP< 3 > . 
1 COORD ( MPO I N , 2 > , LNODS < MELEM , 8 ) ,SHAPE<8) 

GO TO (10,20,30.40) , ISIDE 

SET UP PARAMETERS FOR ELEMENT SIDE WITH FIXED GRADIENTS 

10 CONTINUE 

XEQI V=-l 

YEQI V=P0SGP< IGAUS) 

ITEMP=2 
RTEMP=-1 
GO TO 50 
20 CONTINUE 
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XEQI V=POSGP( IGAUS) 

YEQI V=1 
I TEMP =1 
RTEMP=-1 
GO TO 50 

30 CONTINUE 

XEQI V=1 

YEQI V=POSGP( IGAUS) 

I TEMP =2 
RTEMP=1 
GO TO 50 

40 CONTINUE ' 

XEQ I V=POSGP ( I GAUS ) 

YEQIV=-1 
I TEMP =1 
RTEMP=1 

50 CONTINUE 

DETERMINE SHAPE FUNCTION VALUES AND LOCAL DERIVATIVES 
CALL SHAPES ( DER I V , SHAPE . XEQ I V , YEQI V) 

EVALUATE JACOBIAN AND COSINES OF OUTWARD NORMALS 
IELEMS= IELEM 

CALL D JACOB < COORD , DER I V , DET JB , D JAC I , D JACK , IELEMS , LNODS , 
1 MELEM , MPO I N , NNODP ) 

TEMPY=SQRT ( D JACK ( I TEMP , 1 ) **2 +D JACK ( I TEMP , 2 > **2 ) 
CQSLX=RTEMP*D JACK ( I TEMP , 2 ) /TEMPY 
COSLY = -RTEMP*D JACK ( I TEMP, 1 ) /TEMPY 
SLETH=TEMPY*WE I GP ( IGAUS) 

RETURN 
END 


SUBROUTINE PRESCR ( COORD , EQRHS .GRAOB , IAXSY, IELEM, LNODS .MELEM, 

1 MPO I N , MTOTV . NADFM , NGAUS , NGR AD , NNODP , NODFM . 

2 POSGP , V I SCY , WE I GP ) 

EXTERNAL SUBROUTINES: 

SURF IN 

DIMENSION COORD < MPO I N , 2 ) , EQRHS < MTOTV ) , GRADB ( MPO I N , 4 ) , 

1 LNODS < MELEM , 8 ) . NADFM < MPO IN), NGRAD (MELEM ) . 

2 NODFM (MPO IN) , SHAPE<8 ) , POSGP< 3 ) , WEIGP< 5 ) 

ISIDE=NGRAD( IELEM) 

LOOP TO CARRY OUT GAUSS INTEGRATION 

DO 50 I GAUS =1 .NGAUS 
IELEMS=IELEM 
IGAUSS= IGAUS 

CALL SURF IN(COORD,COSLX, COSLY, IELEMS , I GAUSS , IS IDE , LNODS , 

1 MELEM , MPO I N , NNODP , POSGP , SHAPE , SLETH , WE I GP > 

TDUDN=0 
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RADUS=0 

TDVDN=0 


C 

C 

C 


10 

C 

C 

C 


20 

30 


C 


C 

C 

C 


10 


20 

C 

C 

C 


C 


EVALUATE ( VI SCY*NORMAL-VELOC I TY-GRAD I ENTS ) AT THE GAUSS POINT 

DO 10 INODP=l , NNODP 

KPO I N= I ABS ( LNODS ( I ELEM , I NODP ) ) 

TDUDN=TDUDN+V I SCY*SHAPE ( I NODP)* 

1 ( GRADB ( KPO I N , 1 ) *COSLX+GRADB < KPO I N . 2 ) *COSLY ) 

RADUS =RADUS+SHAPE ( I NODP > *COORD < KPO I N . 2 ) 
TDVDN=TDVDN+VISCY*SHAPE( I NODP)* 

1 ( GRADB ( KPO IN,?) *COSLX+GRADB ( KPO I N , 4 ) *COSLY ) 

CONTINUE 

I F ( I AXSY . EQ . 1 > SLETH=SLETH*RADUS 

INCLUDE BOUNDARY GRADIENT TERMS INTO GLOBAL RHS VECTOR 

DO 20 I NODP =1 .NNODP 

KPO I N= I ABS ( LNODS ( I ELEM, I NODP) ) 

ITOTU=NADFM<KPOIN) 

I TOT V= I TOTU+NODFM < KPO I N ) - 1 

EQRHS ( I TOTU ) =EQRHS ( I TOTU ) +SHAPE ( I NODP ) *TDUDN*SLETH 

EQRHSC ITOTV) = EQRHS < ITOTV)+SHAPE( I NODP ) *TDVDN*SLETH 

CONTINUE 

CONTINUE 

RETURN 

END 

SUBROUT I NE MATR I X ( COORD , DENS Y , EQRHS , FLUMX . I AXSY . LNODS , MELEM . 

1 MPO I N , MTOTV . NADFM , NEVAB , NNODL , NNODP , NODFM . 

2 VARB2 , V I SCY , XFORC , YFORC , VALT , 1 1 TER ) 

DIMENSION AREAW(9 ) , CARLG( 2 , 36 ) , CARPG< 2,72) ,CARTL<2,4) , 

1 CARTP (2.8), COORD ( MPO I N , 2 ) , EQRHS < MTOTV ) , ERHSU ( 8 ) . 

2 ERHS V < 8 ) , FLUMX (20,20), LNODS < MELEM , 8 ) , NADFM ( MPO IN), 

3 NODFM(MPOIN) ,SHALG(?6) ,SHAPG<72) ,SHAPL<4) , 

4 SHAPPC8) .VARB2 (MTOTV) , VALT (MELEM, 9 ) 

INITIALIZE ARRAYS 

DO 10 INODP=l .NNODP 
ERHSU ( I NODP ) = 0 
ERHS V ( I NODP ) = 0 
CONTINUE 

DO 20 IEVAB=1 .NEVAB 
DO 20 JEVAB=1 .NEVAB 
FLUMX ( I E VAB , JEVAB ) = 0 
CONTINUE 

READ IN SHAPE FUNCTION VALUES AND DERIVATIVES 

READ (10) I ELEM , LGAUS . NGAPA , NGAL I , 

1( (CARPG( IDIME, IGAPA) , IDIME=1 ,2) ,SHAPG( IGAPA) , IGAPA=1 , NGAPA) , 

2( ( CARLO ( IDIME, IGALI ) , IDIME=1 .2) ,SHALG( I GAL I ) , IGALI=1 , NGAL I ) , 

3 ( ARE AW ( I GAUS ) , IGAUS=1 , LGAUS) 
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C LOOP TO CARRY OUT GAUSS INTEGRATION 

C 

DO 100 I GAUS =1 , LGAUS 
DAREA=AREAW( IGAUS) 

DO 50 INODP=l .NNODP 

SHAPP< INODP) =SHAPG(NNODP*( I GAUS-1 )+ INODP) 

DO 50 IDIME=1 , 2 

CARTP< 1 0 1 ME , I NODP > =CARPG ( IDIME , NNODP* ( IGAUS-1 )+INODP) 

50 CONTINUE 

DO 40 INODL=l .NNODL 

SHAPL ( I NODL > = SHALG ( NNODL* ( I GAUS - 1 ) + 1 NODL ) 

DO 40 I D I ME = 1,2 

CARTL < I D I ME , I NODL > =CARLG (IDIME. NNODL* ( I GAUS - 1 > + 1 NODL > 

40 CONTINUE 

UVELY=0 
RADUS=0 
VVELY=0 
SRI 1=0 
SR22=0 
SRI 2 = 0 
SREF=0 
VSFAC=1 .0 
C 

C EVALUATE RADIUS AND PREVIOUS VELOCITIES AT GAUSS POINTS 

C 

DO 50 INODP=l .NNOOP 

KPOIN= I ABS(LNODS( IELEM, I NODP) > 

I TOTU =NADFM ( KPO I N ) 

I TOT V = I TOTU+NODFM ( KPO I N ) - 1 
SHAPE =SHAPP< I NODP) 

UVELY=UVELY+VARB2 ( ITOTU)*SHAPE 
RADUS =RADUS +COORD ( KPO I N . 2 ) *SH APE 
VVEL Y = VVEL Y 4-VARB2 < I TOTV ) *SHAPE 
SRI 1 =SR1 1+VARB2( ITOTU)*CARTP( 1 , INODP) 

SR 2 2 = SR 22 +VARB2 < I TOTV ) *C ARTP (2,1 NODP ) 

SR 1 2 = SR 1 2 + . 5* ( VARB2 ( I TOTU > *CARTP (2,1 NODP ) 

1 +VARB2 ( I TOTV ) *CARTP (1,1 NODP ) ) 

50 continue 

SREF=SQRT<2*( SRI 1**2 ) + (SRl 2**2 ) +2*( SR22**2 ) ) 

IF (IITER.EQ.l) GOTO 51 
VSFAC=SREF**- . 95 

51 CONTINUE 

IF( IAXSY.EQ. 1 ) DAREA=DAREA*RADUS 
C 

C PUT BODY FORCES IN LOCAL RHS VECTOR 

C 

DO 60 INODP=l , NNODP 

ERHSU( INODP) =ERHSU( INODP)+SHAPP( INODP)*DAREA*XFORC 
ERHSV( I NODP > =ERHSV ( INODP ) 4-SHAPP( I NODP ) *DAREA*YFORC 
60 CONTINUE 

DO 90 I CONI =1,4 

DO 90 ICON2=l , 2 

IROWU=( I CONI -1 )*5+5*ICON2-2 

I ROWV= I ROWU+ 5 - I CON 2 

INODP=2*( ICON1 -1 )+ICON2 
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SHAP I =SHAPP ( I NODP > 

CARX I =CARTP <1,1 NODP ) 

CARY I =CARTP <2,1 NODP ) 

I ROWP= I ROWU+1 
DO 80 JCON1 =1 , 4 
JCOLP=< JCON1-1 >*5+2 
C 

C PUT PRESSURE TERMS IN MOMENTUM EQUATIONS 

C 

FLUMX< IROWU , JCOLP > =FLUMX( IROWU , JCOLP) +SHAPI *CARTL< 1 . JCQN1 >* 

1 DAREA/DENSY 

FLUMX< IROWV. JCOLP) = FLUMX < IROWV , JCOLP) +SHAPI *CARTL < 2 , JCON1 )* 

1 DAREA/DENSY 

DO 80 JCON2=l , 2 
JCOLU =< JCON1 -2 >*5 + 5* JCON 2-2 
JCOLV= JCOLU+ 3 - JCON2 
JNODP=2*< JCON1-1 )+JCON2 
SHAPJ=SHAPP< JNODP) 

CARXJ=CARTP< 1 , JNODP) 

CARY J =CARTP < 2 , JNODP ) 

C 

C PUT DIFFUSION &. CONVECTION TERMS IN MOMENTUM & ENERGY EQUATIONS 

C 

DIF FU = < CARX I *C ARX J+C AR Y I *CARY J > *V I SCY *DARE A*V ALT < IELEM, IGAUS) 

1 *VSFAC 

CON VC = < UVELY*CARX J+VVELY*CARY J ) *SHAP I *DAREA 
FLUMX < I ROWU , JCOLU ) = FLUMX < I ROWU , JCOLU ) +D I FFU+CONVC 
FLUMX < IROWV , JCOLV ) = FLUMX < I ROWV . JCOLV) +D I FFU+CONVC 
I F < I AXSY . EQ . 1 ) FLUMX < I ROWV , JCOLV ) =FLUMX < I ROWV , JCOLV > + 

1 SHAP I *SHAPJ*VALT< IELEM. IGAUS ) *VI SCY*DAREA/RAOUS**2 

IF < ICON2 . EQ . 2 > GOTO 70 
C 

C FORM CONTINUITY EQUATION 

C 

FLUMX < I ROWP . JCOLU > =FLUMX < I ROWP , JCOLU ) +SH APL < I CON 1 ) *D ARE A*C ARX J 
FLUMX < I ROWP , JCOLV ) =FLUMX < I ROWP , JCOLV ) +SHAPL < I CONI > *DAREA*CARY J 
I F < I AXS Y . EQ . 1 ) FLUMX < I ROWP , JCOLV ) = FLUMX < I ROWP , JCOLV ) + 

1 SHAPL< I CONI >*SHAPJ*DAREA/RADUS 

70 CONTINUE 

80 CONTINUE 

90 CONTINUE 

100 CONTINUE 

C 

C ADD LOCAL RHS VECTOR TO GLOBAL ARRAY 

C 

DO 110 INODP=l ,NNODP 
KPOIN=IABS<LNODS< IELEM, INODP) ) 

I TOTU=NADFM< KPO I N ) 

I TOTV= I TOTU+NODFM < KPO I N > - 1 

EQRHS < I TOTU ) = EQRH S < ITOTU)+ERHSU< INODP) 

EQRHS < I TOTV ) =EQRHS < ITOTV)+ERHSV< INODP) 

110 CONTINUE 

RETURN 
END 
C 
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c 

c 


c 

c 

c 

c 


10 


20 

C 

c 

c 


50 

40 


50 

2000 

1 

2010 


C 

C 

C 


C 

C 

C 

C 


1 

2 

5 

4 

5 


SUBROUT I NE TOLREL < 1 1 TER , MTOTV , NCON V . NTOTV , RELAX . TOLER . VARB 1 . 
VARB2.MPOIN) 

DIMENSION VARB1 (MTOTV) , VARB2(MTOTV> . NADFM(MPOIN) , 
NODFM(MPOIN) 

CHECK TO SEE IF THE SOLUTIONS HAVE CONVERGED 
CHECK FOR CONVERGENCE TO REQUIRED TOLERANCE 

CANLA=0 

NCONV=l 

DO 10 ITOTV=l , NTOTV 

IF(abs<VARBl ( itotv)).LT. .0001 >GO TO 10 

CANGE =ABS( (VARB 1 ( ITOTV) -VARB2 ( I TOTV) ) /VARB1 ( i totv) ) 

I F ( CANGE . GT . TOLER > NCONV= 0 

I F ( CANL A . GT . CANGE > GO TO 10 

CANLA=CANGE 

LTOLA= ITOTV 

CONTINUE 

WRITE (7 , 2000 1LTOLA , CANLA 
IF(NCONV.EQ.O) GO TO 20 
WRITE (7 , 2010 ) 

RETURN 

CONTINUE 

IF( I ITER.EQ. 1 >GO TO 40 

RELAX VARIABLES FOR NEXT ITERATION 

DO 50 I TOT V=1 .NTOTV 

VARB 2 ( I TOTV ) = VARB2 ( I TOTV > +RELAX* ( VARB 1 ( I TOTV > -VARB2 ( I TOTV > ) 

CONTINUE 

RETURN 

CONTINUE 

DO 50 ITOTV=l .NTOTV 
VARB2 ( I TOTV ) = VARB 1 ( I TOTV ) 

CONTINUE 

FORMAT (/57H LARGEST CHANGE OCCURS AT OOF NUMBER . 

14,1 OH CHANGE = .F10.4) 

FORMAT (//45H SOLUTION HAS CONVERGED TO REQUIRED TOLERANCE) 

RETURN 

END 


SUBROUT I NE I TERAT ( BOUDV , COORD , DENSY . EQRHS . GFLUM , GRADB , I AXSY . 
LBOUD . LHEDV , LNODS . MELEM . MFRON , MPO I N . MTOTV . 

NADFM . NBCON , NDOFM , NELEM . NE VAB . NGAUS , NGRAD , 

N I TER , NNODL . NNODP . NODFM , NPO I N . NTOTV . PNORM . 

POSGP , RELAX , TOLER , VARB1 , VARB2 . V I SCY . WE I GP , 

XFORC , YFORC , VALT ) 

EXTERNAL SUBROUTINES: 

PRESCR . FRONTS . WR I TER . TOLREL 
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D I MENS I ON EQRHS < MTOTV > . NGRAD < MELEM ) , COORD ( MPO I N , 2 > . 

1 LNODS ( MELEM , 8 > 

1 , VARB1 (MTOTV) , VARB2.CMTOTV > ,POSGP(3> ,WEIGP<3> , BOUDV ( MTOTV > , 

2 GFLUM < MFRON , MFRON > , LBOUD ( MTOTV ) , LHEDV ( MFRON ) , NODFM ( MPO IN), 

3 PNORM < MFRON ) , NADFM ( MPO IN), GRADB ( MPO I N , 4 ) , VALT ( MELEM , 9 ) 

C 

C SET UP ITERATION COUNTER AND LOOP ADDRESS 

C 

I ITER=0 

10 CONTINUE 

I ITER=I ITER+1 
DO 20 ITOTV=l .NTOTV 
EQRHS <IT0TV)=0 
20 CONTINUE 

C 

C CALL PRESCR FOR ELEMENTS WITH NON-ZERO BOUNDARY GRADIENTS 

C 

DO 30 IELEM=1 ,NELEM 
IF(NGRAD< IELEM) . EQ . 0 )GO TO 30 
IELEMS=IELEM 

CALL PRESCR ( COORD , EQRHS , GRADB , IAXSY , I ELEMS , LNODS . MELEM , 

1 MPO I N , MTOTV , NADFM , NGAUS . NGRAD , NNODP . NODFM , 

2 POSGP , V I SCY , WE I GP ) 

30 CONTINUE 

C 

C 

C CALL FRONTS TO SET UP AND SOLVE GOVERNING EQUATIONS 

C 

CALL FRONTS (BOUDV, COORD, DENS Y, EQRHS. GFLUM. IAXSY. I ITER. 

1 LBOUD , LHEDV , LNODS , MELEM , MFRON , MPO I N , MTOTV , 

2 NADFM , NBCON , NELEM , NEVAB , NNODL , NNODP , NODFM , 

3 NPO IN, NTOTV, PNORM, VARB1 ,VARB2, VI SCY, XFORC, 

4 YFORC, VALT) 
c 

C 

C CALL WRITER TO OUTPUT ITERATION RESULTS 

C 

CALL WR I TER ( 1 1 TER , MPO I N , MTOTV , NADFM , NDOFM , NODFM , NPO I N , 

1 VARB1 , VARB2 ) 

C 

C CALL TOLREL TO CHECK CONVERGENCE. RELAX VALUES IF NOT CONVERGED. 

C 

CALL TOLREL ( 1 1 TER , MTOTV , NCONV , NTOTV , RELAX , TOLER . VARB 1 . VARB2 ) 

C 

C RETURN TO MASTER NUMBER IF ITERATIONS EXCEED MAXIMUM 

C 

I F ( NCONV . EQ . 1 ) RETURN 
IFdITER.LT. NITER) GO TO 10 
WRITE(7 , 2000 ) 

2000 FORMAT(//32H SOLUTION HAS FAILED TO CONVERGE) 

RETURN 

END 

C 

C 

C 
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SUBROUTINE FRONTS ( BOUDV , COORD , DENSY , EQRHS .GFLUM, IAXSY, I ITER, 

1 LBOUD , LHEDV , LNODS , MELEM , MFRON , MPO IN , MTOTV , 

2 NADFM , NBCON , NELEM . NEVAB , NNODL , NNODP , NQDFM , 

? NPO I N , NTOTV , PNORM , VARB1 , VARB2 . V I SCY , XFORC , 

4 YFORC, VALT) 

C 

C EXTERNAL SUBROUTINES: 

C MATRIX 

C 

D I MENS I ON BOUDV ( MTOTV ) , EQRHS < MTOTV ) , FLUMX < 2 0 , 2 0 > , 

1 GFLUM (MFRON, MFRON) .LBOUD (MTOTV) .LHEDV (MFRON) , 

2 LNODS (MELEM, 8 ) , LOCEL( 20 ) , NADFM(MPOIN) ,NDEST(20) , 

3 NODFM( MPO IN), PNORM ( MFRON ) , VARB 1 ( MTOTV ) , VARB2 ( MTOTV > 

4 , COORD ( MPO I N , 2 ) , VALT ( MELEM , 9 ) 

C 

C REWIND TAPES PRIOR TO SOLUTION PROCEDURE 

C 

IF( I ITER.GT. 1 )GO TO 40 
C 

C*** ON FIRST ITERATION ONLY FIND LAST APPEARANCE OF EACH NOD. 

C 

DO 30 IPOIN=l , NPO IN 
LASTE=0 

DO 20 IELEM=1 .NELEM 
DO 10 INODP=l , NNODP 

I F ( LNODS ( I ELEM , I NODP ) . NE . I PO I N > GO TO 10 
LASTE=IELEM 
LASTN= INODP 
GO TO 20 
10 CONTINUE 

20 CONTINUE 

LNODS ( LASTE , LASTN ) = - I PO I N 
30 CONTINUE 

40 CONTINUE 

C 

C*** INITIALISE HEADING AND GRAND FLUID MATRIX. 

C 

REWIND 2? 

REWIND 10 

NCR I T =MFRON-NEVAB 

NFRQN=0 

DO 50 IFRON=l .MFRON 
DO 50 JFRON=l .MFRON 
GFLUM ( I FRON , JFRON ) =0 . 0 
50 CONTINUE 

KELEM=0 
C 

C*** START ASSEMBLY BY FORMING ELEMENTAL MATRICES 

C 

60 CONTINUE 

KELEM=KELEM+ 1 

CALL MATR I X ( COORD , DENSY , EQRHS , FLUMX , I AXSY , LNODS . MELEM , 

1 MPO I N , MTOTV , NADFM .NEVAB , NNODL , NNODP , NODFM , 

2 VARB2 , V I SCY , XFORC , YFORC . VALT , I I TER ) 

KEVAB=0 
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c 

C*** CREAT GLOBAL DOF ARRAY FOR EACH LOCAL ELEMENT DOF. 

C 

DO 70 INODP=l , NNODP 
KPO I N=LNODS ( KELEM , INODP) 

IADFM=NADFM( IABS(KPOIN) > 

LODFM=NODFM( IABS(KPOIN) ) 

DO 70 IODFM=l .LODFM 
KEVAB =KEVAB+1 

LOCEL ( KEVAB ) = I ADFM+ 1 ODFM- 1 
I F < KPO I N . LT . 0 > LOCEL ( KEVAB > = -LOCEL < KEVAB ) 

70 CONTINUE 

C 

C*** FIT EACH DOF INTO THE FRONT WIDTH EXTENDING IF NECESSARY. 

C 

DO 120 IEVAB=1 , NEVAB 
KTOTV=LOCEL( IEVAB) 

I F < NFRON . EQ . 0 ) GO TO 90 
DO 80 IFRON=l .NFRON 
KFRON=IFRON 

I F ( I ABS < KTOTV > . EQ . I ABS ( LHEDV ( KFRON ) > ) GO TO 110 
80 CONTINUE 

90 CONTINUE 

NFRON=NFRON+l 

I F ( NFRON . LE . MFRON > GO TO 100 
WRITE< 7 , 2000) 

STOP 

100 CONTINUE 

NDEST < I EVAB > =NFRQN 
LHEDV (NFRON ) = KTOTV 
GO TO 120 
110 CONTINUE 

NDEST ( I EVAB ) =KFRON 
LHEDV ( KFRON ) =KTOTV 
120 CONTINUE 

C 

C*** ASSEMBLE NEW ELEMENT INTO GRAND FLUID MATRIX. 

C 

DO 150 IEVAB=1 , NEVAB 
IFRON=NDEST ( IEVAB) 

00 150 JEVAB=1 , NEVAB 
JFRON=NDEST < JEVAB ) 

GFLUM( JFRON, IFRON) =GFLUM( OFRQN , I FRON ) +FLUMX ( JEVAB , IEVAB) 

150 CONTINUE 

I F (NFRON.LT. NCR IT. AND. KELEM. LT.NELEM) GO TO 60 
140 CONTINUE 

NFSUM=0 
PI VOT=0 . 0 
C 

C*** CHECK LAST APPEARANCE OF EACH DOF PROCESS BOUNDRY CONDITIONS. 

C 

DO 170 IFRON=l , NFRON 

IF (LHEDV (I FRON) .GE.O)GO TO 170 

NFSUM=1 

IF(LBOUD( IABS(LHEDV( IFRON) ) ) .NE . 1 )GO TO 160 
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260 

270 


280 


290 

500 

C 

C*** 

c 


510 


520 


550 

540 

C 

G*** 

c 


c 

cc*** 

c 


550 


560 


570 

2000 


I F < LP I VT . EQ . 1 > GO TO 270 
DO 260 JFRON=l ,LPIVT-1 

GFLUM< I FRON- 1 , JFRON ) =GFLUM ( I FRON , JFRON ) 

CONTINUE 

CONTINUE 

DO 280 JFRON=LPI VT+1 ,NFRON 


-FACOR*PNORM< JFRON) 


™^i I ^ R0N_1 ’ JFR0N_1 )=GFLUM( I FRON, JFRON) -FACOR*PNORM< JFRON) 
CONTINUE 

ITOTV=IABS<LHEDV< I FRON) ) 

EQRHS< ITOTV)=EQRHS< ITOTV)-FACOR*RHSID 

CONTINUE 

CONTINUE 


WRITE OUT NONFIXED PIVOTAL EQUATION ON TAPE. 

IF <LBOUD<I<TOTV ) .NE.O)GO TO 510 

WR I TE < 2 5 ) NFRON , LP I VT , <LHEDV< I FRON) , PNORM< I FRON) . IFRON=l , NFRON) 
CONTINUE 

DO 520 IFRON=l .NFRON 
GFLUM < I FRON . NFRON >=0.0 
GFLUM < NFRON . I FRON ) = 0 . 0 
CONTINUE 

I F < LP I VT . EQ . NFRON ) GO TO 540 
DO 550 IFRON=LPI VT , NFRON-1 
LHEDV < I FRON ) =LHEDV < I FRON+1 ) 

CONTINUE 
CONTINUE 
NFRON=NFRON- 1 

ASSEMBLE ELIMINATE OR BACK SUBSTITUTE. 

I F < NFRON . GT . NCR I T ) GO TO 140 
I F < KELEM . LT . NELEM > GO TO 60 
I F < NFRON . GT . 0 ) GO TO 140 

BACK SUBSTITUTION 

DO 550 ITOTV=l . NTOTV 
VARB1 < ITOTV)=BOUDV< ITOTV) 

LBOUD < I TOTV > = -LBOUO < I TOTV ) 

CONTINUE 

DO 570 ITOTV=l .NTOTV-NBCON 
BACKSPACE 25 

READ < 2 5 > NFRON . LP I VT , <LHEDV< IFRON) , PNORM< I FRON) . IFRON=l .NFRON) 
KTOTV= I ABS < LHEDV < LP I VT ) ) 

TEMPR=0 . 0 
PNORM<LPIVT)=0. 0 
DO 560 IFRON=l .NFRON 

TEMPR=TEMPR-PNORM< IFR0N)*VARB1 < I ABS < LHEDV < IFRON) ) ) 

CONTINUE 

VARB 1 < KTOTV ) = EQRHS < KTOTV ) +TEMPR 

BACKSPACE 25 

CONTINUE 

RETURN 

FORMAT<//59H PROGRAM HALTED FRONTWIDTH IS TOO SMALL) 
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END 

C 

C 

C 

SUBROUT I NE D INPUT ( BOUDV , COORD , DENS Y , GRADB , I AXSY , LBOUD , LNODS , 
MELEM , MPO I N , MTOTV , NADFM , NBCON , NDOFM , NELEM , 

NE VAB , NGAUS , NGRAD , N I TER , NNODL , NNODP , NODFM , 

NPO I N , NTOTV , RELAX , TOLER . VARB 1 . VARB2 . V I SCY , 

XFORC , YFORC , SLOPE 1 , D , E , FAC , XMAX , YMAX ) 

DIMENSION BOUDV (MTOTV) .COORD (MPO IN, 2) , GRADB < MPO I N , 4 ) , 

1 LBOUD ( MTOTV ) . LNODS ( MELEM , 8 ) , NADFM C MPO IN). 

2 NGRAD (MELEM) .NODFMCMPOIN) ,TITLE(20> , VARB1 (MTOTV) , 

3 VARB2( MTOTV) 

C 

C GIVE VALUES TO THOSE VARIABLES WHICH CANNOT BE CHANGED 

C 

NDOFM=3 
NEVAB=20 
NGAUS =3 
NNODL =4 
NNODP=8 
C 

C READ IN EACH LINE AND ECHO IMMEDIATELY 

C 

READ (40 , 1000)TITLE 
WRITE (7 ,2000 )TITLE 

READ ( 40 , * > I AXSY , N I TER , NEX , NE Y , XMAX , YMAX .angle 
WRITE (7, 2015)1 AXSY, NI TER , NEX , NE Y , XMAX . YMAX .angle 
C READ ( 40 , * ) I AXSY , NELEM , N I TER , NPO I N , NRPON 

C WR I TE ( 7 , 2 Oi 0 > I AXSY , NELEM . N I TER , NPO I N . NRPON 

READ ( 40 , * ) N I CON , NBCON , NEBCN , NNBCN 
WRITEC 7 ,2020 )NICON , NBCON , NEBCN .NNBCN 

CALL MESang ( XMAX , YMAX , NEX , NE Y . LNODS , COORD . NELEM . NPO I N , MELEM 
1 .mpoin, angle) 

C 

C CHECK INITIAL DATA 

C 

CALL D I AGN1 ( NBCON , NEBCN , NELEM , N I CON , NNBCN . NPO I N , NRPON ) 

C 

C READ FLOW PARAMETERS ETC. 

C 

READ ( 40 , * ) DENSY , RELAX , TOLER , V I SCY . XFORC , YFORC 
WR I TE ( 7 , 20 30 > DENSY , RELAX , TOLER , V I SCY , XFORC . YFORC 
READ (40 , *) SLOPE1 , E , D , FAC 
GOTO 25 
C 

C ZERO ALL COORDINATE ARRAY 

C 

DO 10 IPOIN=l , NPO IN 
DO 10 I DIME =1 , 2 
COORD ( IPOIN, IDIME)=0 
10 CONTINUE 

C 

C READ IN NODAL COORDINATES (X THEN Y) 

C 
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DO 20 IPOIN=l .NRPON 

READ( 40 , *) JPOIN, (COORD( JPOIN , I DIME) , IDIME = 1 ,2) 

20 CONTINUE 

C 

C READ IN NODAL CONNECTIONS AND ECHO IMMEDIATELY 

C 

25 CONTINUE 

WRITEC7 , 2040) 

DO 50 IELEM=1 .NELEM 
NGRAD ( I ELEM > = 0 

C READ (40,*) I ELEM , ( LNODS ( I ELEM , I NODP ) , I NODP= 1 , NNODP ) 

WRITE <7, 2050) I ELEM, (LNODS ( I ELEM , INODP) , INODP=l , NNODP) 

50 CONTINUE 

C 

C SET UP VECTOR GIVING THE D.O.F. AT EACH NODE 

C 

ITEMP=NNODP-l 
DO 40 IELEM=1 , NELEM 
DO 40 I NODP =1 . I TEMP, 2 
NODFM ( LNODS ( I ELEM , I NODP ) ) =NDOFM 
NODFM ( LNODS ( I ELEM, INOOP+1 ) )=NDOFM-l 
40 CONTINUE 

GOTO 71 
C 

C INTERPOLATE MIDSIDE NOOE COORDINATES WHERE NOT GIVEN 

C 

DO 70 IELEM=1 .NELEM 
DO 60 INODP=2 , NNODP , 2 
I PO I N=LNODS ( I ELEM . I NODP ) 

TEMPY= ABS ( COORD ( I PO I N , 1 ) ) +ABS ( COORD ( I PO I N . 2 ) ) 

IF(TEMPY.NE.O.O) GO TO 60 
JPOIN=LNODS( I ELEM, INODP-1 > 

KNODP= INOOP+1 
I F ( KNODP . GT . NNODP ) KNODP = 1 
KPO I N=LNODS ( I ELEM , KNODP ) 

DO 50 IDIME=1 , 2 

COORD ( I PO IN, I D I ME > = ( COORD ( JPO I N , I D I ME ) +COORD ( KPO I N , IDIME) )*. 5 
50 CONTINUE 

60 CONTINUE 

70 CONTINUE 

71 CONTINUE 
C 

C SET UP VECTOR WITH THE FIRST D.O.F. AT EACH NODE 

C 

NADFM( 1 ) =1 

DO 80 IPOIN=2 .NPOIN 

NADFM( IPO IN) =NADFM( IPOIN-1 )+NODFM( IPOIN-1 ) 

80 CONTINUE 

NTOTV=NADFM ( NPO I N ) +NODFM ( NPO I N ) - 1 
C 

C INITIALIZE REMAINING ARRAYS 

C 

DO 90 ITOTV=l ,NTOTV 
LBOUD ( I TOT V > = 0 
BOUDV ( I TOT V > = 0 
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VARB2( ITOTV)=0 
90 CONTINUE 

C 

C READ IN INITIAL CONDITIONS (GUESS VALUES) IF ANY 

C 

I F (NICON . EQ . 0 > GO TO 110 
DO 100 I ICON=l , NICON 
READ (40 ,*) IPOIN , IDOFM.TEMPY 
I F ( I DOFM . GT . 1 > I DOFM=NODFM (IPOIN)-l 
JTOTV=NADFM ( I PO I N ) + 1 DOFM- 1 
VARB1 ( JTOTV)=TEMPY 
VARB2 ( JTOTV) =TEMPY 
100 CONTINUE 

110 CONTINUE 

C 

C WRITE OUT COORDINATES AND ARRAYS SHOWING DOF 

C 

WRITE (7 ,2060) 

DO 130 IPOIN=l ,NPOIN 

IF(NODFM( IPOIN) .NE.NDOFM) GO TO 120 

WRITE(7. 2070) IPOIN, NADFM( IPOIN) ,NODFM( IPOIN) , 

1 ( COORD ( IPOIN, IDIME) , IDIME = 1 ,2) , 

2 ( VARB1 (NADFM( I POIN) -1 + IODFM) , IODFM=l ,NODFM( IPOIN) ) 

GO TO 130 

120 CONTINUE 

WRITE(7 , 2080) I POIN, NADFM( IPOIN) , NODFM( IPOIN) . 

1 ( COORD ( IPOIN, IDIME) , IDIME=1 ,2) , 

2 (VARB1 (NADFM( IPOIN) -1 + IODFM) , IODFM=l ,NODFM( IPOIN) ) 

130 CONTINUE 

C 

C CHECK NODAL CONNECTIONS AND COORDINATES 

C 

CALL D I AGN2 ( COORD , DENS Y , LNODS , MELEM , MPO I N , NELEM , N I CON , 
1 NNODP , NPO I N , NTOTV , V I SCY ) 

C 

C READ IN BOUNDARY CONDITIONS 

C 

WRITE (7 ,2090) 

DO 170 IBCON=l .NBCON 
IF (NBCON. EQ.O) GOTO 171 
READ (40 ,*> IPOIN , IDOFM, BVALU 
ITOTV=NADFM(IPOIN)+IDOFM-l 
I F ( I DOFM . GT . NODFM ( I PO I N ) ) I TOTV= I TOT V- 1 
LBOUD ( I TOTV > = 1 
BOUD V ( I TOTV ) = BVALU 
GO TO (140,150,160) , IDOFM 
140 CONTINUE 

WRITE(7, 2100) IPOIN, BVALU 
GO TO 170 
150 CONTINUE 

WR I TE (7, 2110)1 POIN, BVALU 
IF(NODFM( IPOIN) .NE. IDOFM) GO TO 170 
WR I TE (7,2120) I BCON , IPOIN 
STOP 
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160 

170 

171 
C 


C 

C 


180 


190 

C 

C 

C 


200 


210 


220 


2 30 

240 

1000 

2000 

2010 


2015 


2020 

2050 


2040 


CONTINUE 

WRITE <7 , 2150)1 PO IN, BVALU 

CONTINUE 

CONTINUE 


READ IN ELEMENTS AND SIDES WITH GRADIENT BOUNDARY CONDITIONS 

IF(NEBCN.EQ.O) RETURN 
WRITE <7 ,2140) 

00 180 IEBCN=1 .NEBCN 

READ < 4 0 . * ) I ELEM . NGRAD < I ELEM ) 

WRITE (7 ,2150) I ELEM, NGRAD ( I ELEM) 

CONTINUE 
WRITE (7, 2160) 

DO 190 IPOIN=l .NPOIN 
DO 190 IGRAD=1 , 4 
GRADB( IPO IN, IGRAD) =0 
CONTINUE 

READ IN VALUES OF THE GRADIENTS ON BOUNDARY NODES 
DO 240 INBCN=1 .NNBCN 

READ ( 40 , *• ) I PO I N , I GRAD , GRADB ( I PO I N , I GRAD) 

GO TO (200,210,220,250), IGRAD 
CONTINUE 

WR ITE(7, 2100)1 PO IN, GRADB ( I PO I N , I GRAD ) 

GO TO 240 
CONTINUE 

WRITE(7 . 2110) IPOIN.GRADB< IPOIN , IGRAD) 

GO TO 240 
CONTINUE 

WRITE< 7. 21 50) IPOIN, GRADB (IPOIN. IGRAD) 

GO TO 240 
CONTINUE 

WRITE (7, 21 70) IPOIN, GRADB (IPOIN, IGRAD) 

CONTINUE 
RETURN 
FORMAT (80A1 ) 

FORMAT ( 1 4H OUTPUT FILE :,//lX,20A4) 

FORMAT (//l 5H CONTROL DATA, /I 5H **#***♦*****,// 

18H IAXSY = , 1 4 , 4X , 8H NELEM =,I4,4X,8H NITER =,I4,4X, 

28H NPOIN = , 1 4 , 4X , 8H NRPON =,I4) 

FORMAT (//l 5H CONTROL DATA , /l 5H ************,// 

18H IAXSY = , 1 4 , 6X , 8H NITER =,I4,6X,8H NEX =,I4,6X, 

28H NEY =,I4,/8H XMAX =,F6.2,4X,8H YMAX =,F6.2,4X, 

58H angle =,f7.2) 

FORMAT ( 8H NICON =.I4,6X,8H NBCON =, 

1I4.6X.8H NEBCN =,I4,6X,8H NNBCN =,I4) 

FORMAT ( / /20H PHYSICAL PROPERTI ES , /20H *******************,// 
18H DENSY = , F 1 0 . 5 , 4X , 8H RELAX = , FI 0 . 5 , 4X , 8H TOLER =, 

2F10.5./8H VISCY = , FI 0 . 5 , 4X , 8H XFORC =,F10.5.4X. 

58H YFORC =,F10.5) 

FORMAT ( / / 1 7H ELEMENT TOPOLOGY , /I 7H ****************,// 
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18H ELEMENT. 6X. 12HN0DE NUMBERS) 

2050 FORMATdX. I5.6X.8I5) 

2060 F0RMAT(//16H NODE POINT DATA./16H *******##*##***,// 

134H NODE NADFM NODFM X-COORD Y-C00RD.5X, 

225HU-VEL PRESSURE V-VEL) 

2070 FORMAT (3<I4,1X),5(F9.4,1X)> 

2080 F0RMAT<5( 14, IX) , 3 (F9 . 4 , IX) , 1 OX , F9 . 4 ) 

2090 FORMAT(//20H BOUNDARY CONDITIONS , /20H #*##*•*#***##*###*#* , / / 
135H NODE U-FIXED PRESSURE V-FIXED ) 

2100 FORMAT < 1 5 , F 9 . 4 ) 

2110 FORMAT < 1 5 , 1 OX , F9 . 4 > 

2120 FORMAT < /27H ERROR!! BOUNDARY CONDITION. 14 . 

142H DEFINES A PRESSURE AT MIDSIDE NODE NUMBER, 14) 

2130 FORMAT <15, 2 OX, F9. 4) 

2140 FORMAT < /25H VELOCITY GRADIENT VALUES ,/. IX , 24< '*') . 

1//14H ELEMENT SIDE) 

2150 FORMAT < 1 5 , 3X , 1 5 ) 

2160 FORMAT ( // 44H NODE DU/DX DU/DY DV/DX DV/DY) 

2170 FORMAT <15, 30X ,F9 - 4) 

END 


SUBROUT I NE D I AGN 1 < NBCON , NEBCN , NELEM , N I CON , NNBCN , NPO I N , NRPON > 
DIMENSION NECHO (80), NEROR < 8 ) 

SCRUTINY OF CONTROL DATA 

DO 10 IEROR=l , 8 
0 NEROR < I EROR > =0 

SCRUTINIZE CONTROL DATA AND PRINT ERROR MESSAGE 

IF(NPOIN.LT.O) NEROR ( 1 >=1 

COMMENT NEXT LINE IF NOT READING NRPON. 

I F ( NRPON . LE . 0 . OR . NRPON . GT . NPO I N ) NEROR ( 2 ) = 1 
IF(NELEM.LE.O) NEROR(3)=l 

I F ( NPO I N . GT . NELEM*8 ) NER0R<4)=1 

I F < NBCON . GT . NPO I N*3 ) NEROR(5) = l 

IF(NICON.GT.NPOIN*3> NEROR(6)=l 

IF<NEBCN.GT. NELEM) NER0R<7)=1 

I F < NNBCN . GT . NPO I N*4 ) NER0R<8) = 1 

CHECK ON ERROR. IF UNITY, PRINT. 

JEROR=0 

DO 20 IEROR=l , 8 

I F < NEROR < I EROR ) . EQ . 0 ) GO TO 20 

JEROR=l 

WRITE OUT ERROR NUMBER 

WR I TE <7, 2000) I EROR 
20 CONTINUE 

IF ( JEROR. EQ. 0) RETURN 
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c 

C LIST DATA REMAINING AFTER CONTROL PARAMETERS 

C 

WRITE(7, 2010> 

50 CONTINUE 

READ(40 , 1 000 >NECHO 
WRITE (7, 2020) NECHO 
GO TO 50 

1000 FORMAT < 80A1 ) 

2000 FORMAT ( / 1 1 OX , 5 5HCONTROL PARAMETER ERROR*****ERROR .15) 

2010 FORMAT(/10X. 57HDATA FOLLOWING ERROR IN CONTROL CARDS/) 

2020 FORMAT ( 1 OX , 80A1 ) 

END 

C 

C 

C 

SUBROUT I NE D I AGN2 ( COORD . DENSY . LNODS , MELEM , MPO I N , NELEM , N I CON , 
1 NNODP , NPO I N , NTOTV , V I SCY ) 

DIMENSION COORD(MPOIN, 2 ) , LNODS (MELEM, 8 ) .NECHO (80) ,NEROR<15) 

C 

C SCRUTINIZE ELEMENT AND NODAL POINT DATA 

C INITIALIZE ERROR ARRAY 

C 

DO 10 I EROR= 1 ,15 
NEROR < I EROR ) = 0 
10 CONTINUE - 

C 

C CHECK PHYSICAL PROPERTIES 

C 

IF(DENSY.LE.O .OR.VISCY.LE.O) NER0R<9)=1 
C 

C CHECK NODAL COORDINATES 

C CHECK FOR TWO IDENTICAL COORDINATES 

C 

00 40 IPOIN=2 , NPO IN 

JPOIN=IPOIN-l 

DO 50 KPOIN=l , JPOIN 

I F < COORD < I PO I N , 1 > . NE . COORD C KPO I N , 1 ) .OR.COORD( IPOIN.2) . 

1 NE . COORO ( KPO I N , 2 > ) GO TO 20 



NEROR ( 1 0 > = 1 


20 

CONTINUE 


50 

CONTINUE 


40 

r 

CONTINUE 


u 

c 

CHECK ELEMENT NODAL 

NUMBERS 

c 

A - REPETITION OF A 

NODE NUMBER IN ONE ELEMENT 

c 

r- 

B - NODE NUMBER OUTSIDE PERMISSIBLE BOUNDS 


DO 50 IPOIN=l , NPOIN 
DO 50 IELEM=1 , NELEM 
LCONT =0 

DO 50 I NODE =1 .NNODP 



I F ( LNODS ( I ELEM , I NODE ) . NE . I PO I N ) GO TO 50 
LCONT =LCONT + 1 

IF(LCONT.GT.l) NEROR(ll)=l 
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30 CONTINUE 

DO 60 IELEM=1 .NELEM 
DO 60 I NODE =1 .NNODP 

I F < LNODS ( I ELEM , INODE) .LE . 0 . OR .LNODS( IELEM, INODE) .GT. NPOIN) 
1NER0R( 1 2 > = 1 
60 CONTINUE 

C 

C CHECK ON NUMBER OF INITIAL CONDITIONS INPUT 

C 

I F < N I CON . GT . NTOTV ) NER0R<13>=1 
C 

C WRITE OUT ERROR NUMBER 

C 

JEROR=0 

DO 70 IEROR=9 , 1 3 

I F < NEROR < I EROR ) . EQ . 0 ) GO TO 70 

JEROR=l 

WRITE<7 , 2000 ) I EROR 
70 CONTINUE 

C 

C DECIDE WHETHER TO CONTINUE OR ECHO DATA 

C 

IF(JEROR.EQ.O) RETURN 
C 

C LIST REMAINING DATA 

C 

WRITE (7,2010) 

80 CONTINUE 

READ( 40 , 1 000 ) NECHO 
WRITE(7, 2020) NECHO 
GO TO 80 

1000 FORMAT(80A1 ) 

2000 FORMAT ( / 1 OX , 1 OHDATA ERROR ,15) 

2010 FORMAT ( / 1 OX , 1 4HREMA I N I NG DATA/) 

2020 FORMAT (1 OX, 80 Al) 

END 

C 

SUBROUT I NE MATPRT < COORD , VALT , LNODS . MELEM , MPO I N , NELEM , NGAUS . 
1 POSGP , XMAX , YMAX , SLOPE 1 , D , E , F AC , NNODP ) 

D I MENS I ON COORD < MPO I N . 2 > , VALT < MELEM , 9 ) , POSGP ( 3 ) , 

1 SHAPE<8> , LNODS <MELEM, 8 ) .DERI V< 2 , 8 ) 

LGAU S = NGAUS *NGAU S 
DO 11 IELEM=1 .NELEM 
DO 10 KGAUS=1 .LGAUS 
VALT < I ELEM , KGAUS > = 1.0 

10 CONTINUE 

11 CONTINUE 

DO 100 IELEM=1 .NELEM 
KGAUS =0 

DO 90 IGAUS=1 .NGAUS 
DO 90 JGAUS=1 .NGAUS 
KGAUS =KGAUS+1 

C FINDING COORDINATES OF GAUSS-POINTS 

XEQI V=POSGP< IGAUS) 
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CALL SHAPE 8 < DER I V , SHAPE , XEQ I V . YEQ I V ) 

XC=0.0 

YC=0.0 

c shape-f unc=l , 8 

DO 80 I NODP =1 , NNODP 

KPO I N= I ABS ( LNODS < I ELEM , I NODP > > 

XC=XC+SHAPE ( I NODP ) *COORD < KPO I N , 1 ) 

YC= YC+SHAPE < I NODP > *COORD ( KPO I N , 2 > 

80 CONTINUE 

C CHECK THE POSITION OF GAUSS-POINTS 

YCRU=SLOPE 1 * ( XC-XMAX/ 2 ) + ( YMAX/2 +D*SQRT ( 1 +SLOPE1 **2 > > 
YCRL=SL0PE1 *( XC-XMAX/2 > + ( YMAX/2 -D*SQRT( 1 +SLOPE1 **2 ) > 
SLOPE2=-l /SLOPE1 

XCRU= ( YC-YMAX/2 > / SLOPE 2 + ( XMAX/ 2 +E*SQRT ( 1 +SLOPE 1 **2 ) > 

XCRL= (YC- YMAX/2 ) /SLOPE2+( XMAX/2-E*SQRT ( 1 +SLOPE1 **2 > > 

IF ((XC.LE.XCRU) . AND . < XC . GE . XCRL > . AND . ( YC . LE . YCRU) . 

1 AND. (YC.GE.YCRL) ) THEN 
VALT ( I ELEM , KGAUS > = FAC 
ELSE 
END IF 

90 CONTINUE 

100 CONTINUE 

DO 110 IELEM=1 .NELEM 

WRITE<7,*> I ELEM, ( VALT ( I ELEM, KGAUS) , KGAUS = 1 , LGAUS ) 

110 CONTINUE 

RETURN 
END 

C 

C 

SUBROUT I NE WR I TER < 1 1 TER , MPO I N , MTOTV , NADFM . NDOFM . NODFM , NPO I N , 

1 VARB1 , VARB2) 

DIMENSION NADFM ( MPO IN), NODFM < MPO IN) .VARB1 (MTOTV) , 

2 VARB2 (MTOTV ) 

WRITE(7 , 2000)1 ITER 
DO 20 IPOIN=l ,NPOIN 
IODFM=NODFM( IPOIN) 

IADFM=NADFM( IPOIN) 

IFdODFM.EQ. NDOFM) GO TO 1 0 

c WRITE(7, 2010) IPOIN, (VARB1( I ADFM+ JODFM-1 ) . JODFM=l . IODFM) , 

c 1 ( VARB2 ( I ADFM+ JODFM- 1 ) . JODFM= 1 , I ODFM > 

WRITE(7 , 2010 ) IPOIN , (VARB1 ( I ADFM+ JODFM-1 ) , JOOFM=l , IODFM) 

GO TO 20 
10 CONTINUE 

C WRITE(7, 2020) IPOIN, (VARB1( IADFM+ JODFM-1 ) , JODFM=l , IODFM) , 

C 1 (VARB2( I ADFM+ JODFM-1 ), JODFM=l . IODFM) 

WRITE(7, 2020) IPOIN, (VARB1 ( I ADFM+ JODFM- 1 ) . JODFM=l , IODFM) 

20 CONTINUE 

RETURN 

2000 FORMAT ( //29H RESULTS FOR ITERATION NUMBER, 14. 

1 //14X , 5 ( 5HNEW , 7X) , 4X , 3 ( 5HOLD , 7X > / , 

27H NODE , 2(4X , 30HU-VELOC ITY PRESSURE V-VELOCITY)) 

2010 FORMAT (16, 4X , ( E 1 0 . 3 . 10X.E10. 5, 5X) ) 

2020 FORMAT (16, 4X , ( 5E1 0 . 5 , 5X) ) 
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c 

c 

SUBROUT I NE MESang ( XMAX , YMAX , NEX , NE Y , LNODS , COORD , NELEM , NPO I N , 
1 MELEM.MPOIN, angle) 

D I MENS I ON LNODS < MELEM , 8 > , COORD ( MPO I N , 2 > , b ( mpo in) 
DX=XMAX/(2*NEX) 

DY=YMAX/(2*NEY) 

NPO I N= < 2*NEX+ 1 ) * ( 2*NEY+1 ) -NEX*NEY 

NELEM=NEX*NEY 

a=ymax/tand<angl e ) 

i f < 2*a . 1 t . xmax) goto 5 

wr i te( 5 , 1000 > 

1000 f ormat < //10x , ' i ncompat i bl e mesh geometry’) 
stop 

5 dxt=(xmax-2*a)/(2*nex> 

ddx=dx-dxt 
do 25 ipoin=l .npoin 
b( ipoin)=0 
25 continue 

DO 50 IELEM=1 , NELEM 
J1 = I ELEM/NEX 

IF< IELEM. EQ. J1*NEX) GOTO 10 

J=J1+1 

GOTO 20 

10 CONTINUE 

J=J1 

20 CONTINUE 

I = IELEM- ( J-l )*NEX 
N1=<3*NEX+2)*(J-1)+2*<I-1)+1 
LNODS ( IELEM, 1>=N1 
LNODS < IELEM. 2>=N1+1 
LNODS < I ELEM , 5 > =N1 +2 
LNODS < I ELEM , 4 ) =2*NEX- I +N1 +3 
LNODS ( I ELEM , 5 ) = 3*NEX+N1 +4 
LNODS < I ELEM , 6 > =LNODS ( I ELEM , 5 ) - 1 
LNODS ( I ELEM , 7 ) =LNODS < I ELEM , 6 ) - 1 
LNODS < I ELEM , 8 ) =LNODS ( I ELEM , 4 ) - 1 
X1=(I-1)*2*DX 
Yl = < JJ-l )*2*DY 
COORD < LNODS ( I ELEM . 1 ) . 1 ) =X1 
COORD < LNODS < I ELEM , 2 > . 1 > =X1 +DX 
COORD ( LNODS ( IELEM, 3 ) , 1 ) =X1 +2*DX 
COORD < LNODS < I ELEM , 4 ) , 1 ) = X 1 + 2 *DX 
COORD ( LNODS <IELEM,5),l)=Xl+2 *DX 
COORD < LNODS < IELEM, 6) . 1 )=X1+DX 
COORD < LNODS (IELEM,7),1)=X1 
COORD < LNODS < IELEM, 8) , 1 )=X1 
COORD < LNODS < IELEM, 1 ) ,2)=Y1 
COORD(LNODS< IELEM , 2 > , 2 ) =Y1 
COORD < LNODS < I ELEM , 3 ) , 2 ) = Y1 
COORD (LNODS ( IELEM , 4) , 2 ) =Y1 +DY 
COORD ( LNODS ( I ELEM , 5 ) , 2 ) =Y1 +2*DY 
COORD(LNODS( I ELEM, 6) , 2)=Y1+2*DY 
COORD ( LNODS ( I ELEM , 7 ) , 2 ) =Y1 +2*DY 
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tanl=ymax/<a-2*( i -1 >*ddx> 
tan2=ymax/(a-<2*i -1 >*ddx> 
tan3=ymax/ (a-2*i *ddx> 

b< 1 nods < i e 1 em, 1 > >= coord < lnods(ielem,l),2)/tanl 
b( lnods( i e 1 em, 2 > > = coord ( lnods ( i e lem, 2> , 2> /tan2 
b< lnods< i e 1 em, 3 > >=coord< lnods < i elem, 3) , 2 >/tan3 
b< lnods < i e 1 em, 4 > )= coord ( lnods ( ielem,4>,2>/ tan 3 
b< lnods < i e 1 em, 5 > > = coord < lnods ( i e 1 em, 3 > , 2 > /tan 3 
b< lnods < i e 1 em, 6 > > = coord ( lnods (ielem,6>,2)/tan2 
b< lnods < i e 1 em, 7> ) = coord < lnods (ielem,7>,2)/tanl 
b( 1 nods <ielem,8))=coord( lnods (ielem,8),2>/tanl 
CONTINUE 

do 40 i po i n=l , npo i n 

coord ( i po i n , 1 > =coord< i po i n , 1 ) +b( i po in) 

cont i nue 

RETURN 

END 
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OUTPUT FILE 


a 

0 

5 

V 

f 

4 

CONTROL 

DATA 





************ 





IAXSY = 

0 

NITER 

s 

4 

NEX 

XMAX = 

1 .00 

YMAX 

= 

1.00 

angle 

NICON = 

0 

NBCON 

= 

86 

NEBCN 


PHYSICAL PROPERTIES 
******************* 


DENSY = 1.00000 RELAX 

VISCY = 100.00000 XFORC 


.50000 TOLER = 

.00000 YFORC = 


ELEMENT TOPOLOGY 
**************** 


ELEMENT 

NODE 

NUMBERS 






l 

l 

2 

5 

19 

29 

28 

27 

18 

2 

5 

4 

5 

20 

51 

50 

29 

19 

5 

5 

6 

7 

21 

55 

52 

51 

20 

4 

7 

8 

9 

22 

55 

54 

55 

21 

5 

9 

10 

11 

25 

57 

56 

55 

22 

6 

11 

12 

15 

24 

59 

58 

57 

25 

7 

15 

14 

15 

25 

41 

40 

59 

24 

8 

15 

16 

17 

26 

45 

42 

41 

25 

9 

27 

28 

29 

45 

55 

54 

55 

44 

10 

29 

50 

51 

46 

57 

56 

55 

45 

11 

51 

52 

55 

47 

59 

58 

57 

46 

12 

55 

54 

55 

48 

61 

60 

59 

47 

1? 

55 

56 

57 

49 

65 

62 

61 

48 

14 

57 

58 

59 

50 

65 

64 

65 

49 

15 

59 

40 

41 

51 

67 

66 

65 

50 

16 

41 

42 

45 

52 

69 

68 

67 

51 

17 

55 

54 

55 

71 

81 

80 

79 

70 

18 

55 

56 

57 

72 

85 

82 

81 

71 

19 

57 

58 

59 

75 

85 

84 

85 

72 

20 

59 

60 

61 

74 

87 

86 

85 

75 

21 

61 

62 

65 

75 

89 

88 

87 

74 

22 

65 

64 

65 

76 

91 

90 

89 

75 

25 

65 

66 

67 

77 

95 

92 

91 

76 

24 

67 

68 

69 

78 

95 

94 

95 

77 

25 

79 

80 

81 

97 

107 

106 

105 

96 

26 

81 

82 

85 

98 

109 

108 

107 

97 

27 

85 

84 

85 

99 

111 

110 

109 

98 

28 

85 

86 

87 

100 

115 

112 

111 

99 

29 

87 

88 

89 

101 

115 

114 

115 

100 

50 

89 

90 

91 

102 

117 

116 

115 

101 

51 

91 

92 

95 

105 

119 

118 

117 

102 

52 

95 

94 

95 

104 

121 

120 

119 

105 

55 

105 

106 

107 

125 

155 

152 

151 

122 

54 

107 

108 

109 

124 

155 

154 

155 

125 

55 

109 

110 

111 

125 

157 

156 

155 

124 

56 

111 

112 

115 

126 

159 

158 

157 

125 


57 

115 

114 

115 

127 

141 

140 

159 

126 

58 

115 

116 

117 

128 

145 

142 

141 

127 

59 

117 

118 

119 

129 

145 

144 

143 

128 

40 

119 

120 

121 

150 

147 

146 

145 

129 


1 / d = 1 . 5 

NEY = 8 

NNBCN = 0 

.01000 
. 00000 



41 

1 >1 

132 

1 33 

14 ? 

17 ? 

1 >« 

IP / 

i 4 b 

42 

133 

134 

133 

130 

161 

160 

139 

149 

43 

133 

136 

137 

131 

163 

162 

161 

130 

44 

137 

138 

139 

132 

163 

164 

163 

131 

43 

139 

140 

141 

133 

167 

166 

163 

132 

46 

141 

142 

143 

134 

169 

168 

167 

133 

47 

143 

144 

143 

133 

171 

170 

169 

134 

48 

143 

146 

147 

136 

173 

172 

171 

133 

49 

137 

138 

139 

173 

183 

184 

183 

174 

30 

139 

160 

161 

176 

187 

186 

183 

173 

31 

161 

162 

163 

177 

189 

188 

187 

176 

32 

163 

164 

163 

178 

191 

190 

189 

177 

33 

163 

166 

167 

179 

193 

192 

191 

178 

34 

167 

168 

169 

180 

193 

194 

193 

179 

33 

169 

170 

171 

181 

197 

196 

193 

180 

36 

171 

172 

173 

182 

199 

198 

197 

181 

37 

183 

184 

183 

201 

211 

210 

209 

200 

38 

183 

186 

187 

202 

213 

212 

211 

201 

39 

187 

188 

189 

203 

213 

214 

213 

202 

60 

189 

190 

191 

204 

217 

216 

213 

203 

61 

191 

192 

193 

203 

219 

218 

217 

204 

62 

193 

194 

193 

206 

221 

220 

219 

205 

63 

193 

196 

197 

207 

223 

222 

221 

206 

64 

197 

198 

199 

208 

223 

224 

223 

207 


NODE 

POINT 

DATA 






*** *********** * 






NODE 

NADFM 

NODFM 

X-COORD 

Y-COORD 

U-VEL 

p 

V-VEL 

1 

1 

5 

. 0000 

.0000 

.0000 

.0000 

.0000 

2 

4 

2 

.0625 

.0000 

.0000 


.0000 

3 

6 

3 

. 1250 

.0000 

.0000 

.0000 

.0000 

4 

9 

2 

.1875 

.0000 

.0000 


.0000 

5 

11 

3 

.2300 

.0000 

.0000 

.0000 

.0000 

6 

14 

2 

.3125 

.0000 

.0000 


.0000 

7 

16 

3 

. 3750 

.0000 

.0000 

.0000 

.0000 

8 

19 

2 

.4375 

.0000 

.0000 


.0000 

9 

21 

3 

. 5000 

.0000 

.0000 

.0000 

.0000 

10 

24 

2 

.5625 

.0000 

.0000 


.0000 

11 

26 

3 

.6250 

.0000 

.0000 

.0000 

.0000 

12 

29 

2 

.6875 

.0000 

.0000 


.0000 

13 

31 

3 

.7500 

.0000 

.0000 

.0000 

.0000 

14 

34 

2 

.8125 

.0000 

.0000 


.0000 

13 

36 

3 

.8750 

.0000 

.0000 

.0000 

.0000 

16 

39 

2 

.9375 

.0000 

.0000 


.0000 

17 

41 

3 

1.0000 

.0000 

.0000 

.0000 

.0000 

18 

44 

2 

.0000 

.0625 

.0000 


.0000 

19 

46 

2 

.1250 

.0625 

.0000 


.0000 

20 

48 

2 

.2500 

.0625 

.0000 


.0000 

21 

30 

2 

.3730 

.0625 

.0000 


.0000 

22 

32 

2 

.3000 

.0625 

.0000 


.0000 

25 

54 

2 

.6250 

.0625 

.0000 


.0000 

24 

56 

2 

.7500 

.0625 

.0000 


.0000 

23 

38 

2 

.8730 

.0625 

.0000 


.0000 

26 

60 

2 

1.0000 

.0625 

.0000 


.0000 

27 

62 

3 

.0000 

.1250 

.0000 

.0000 

.0000 

28 

63 

2 

.0625 

.1250 

.0000 


.0000 

29 

67 

3 

.1230 

.1250 

.0000 

.0000 

.0000 

30 

70 

2 

.1875 

.1250 

.0000 


.0000 

51 

72 

3 

.2500 

.1250 

.0000 

.0000 

.0000 

32 

75 

2 

.3125 

.1250 

.0000 


.0000 


33 

77 

3 

.3750 

. 1250 

.0000 

.0000 

. 0000 

34 

80 

o 

i. 

.4375 

.1250 

.0000 


.0000 

35 

82 

3 

.5000 

.1250 

.0000 

.0000 

.0000 

36 

8 5 

2 

.5625 

.1250 

.0000 


. 0000 



37 8 / > 

58 90 2 

59 92 J 

40 9 $ 2 

41 97 j 

42 100 2 

43 102 3 

44 103 2 

43 107 2 

46 109 2 

47 111 2 

48 113 2 

49 113 2 


50 

117 

2 

51 

119 

2 

52 

121 

2 

53 

123 

3 

54 

126 

2 

55 

128 

3 

56 

131 

2 

57 

133 

3 

58 

136 

2 

59 

138 

3 

60 

141 

2 

61 

143 

3 

62 

146 

2 

63 

148 

3 

64 

151 

2 

65 

153 

3 

66 

156 

2 

67 

158 

3 

68 

161 

2 

69 

163 

3 

70 

166 

2 

71 

168 

2 

72 

170 

2 

73 

172 

2 

74 

174 

2 

75 

176 

2 

76 

178 

2 

77 

180 

2 

78 

182 

2 

79 

184 

3 

80 

187 

2 

81 

189 

3 

82 

192 

2 

83 

194 

3 

84 

197 

2 

85 

199 

3 

86 

202 

2 

87 

204 

3 

88 

207 

2 

89 

209 

3 

90 

212 

2 

91 

214 

3 

92 

217 

2 

93 

219 

5 

94 

222 

2 

95 

224 

3 

96 

227 

2 

97 

229 

2 

98 

231 

2 


. 62 >U 

.6873 
.7300 
.8123 
.8730 
.9373 
1.0000 
.0000 
.1230 
.2300 
. 3730 
.3000 
.6230 
.7300 
.8730 
1 .0000 
.0000 
.0623 
.1230 
.1873 
.2300 
. 3123 
.3730 
.4373 
.3000 
.3623 
.6230 
.6873 
.7300 
.8123 
.8730 
.9373 
1.0000 
.0000 
.1230 
.2300 
.3730 
.3000 
.6230 
.7300 
: 8730 
1 .0000 
.0000 
.0623 
.1230 
.1873 
.2300 
.3123 
.3730 
.4373 
.3000 
.3623 
.6230 
.6873 
.7300 
.8123 
.8730 
.9373 
1 .0000 
.0000 
.1230 
.2300 


.1230 
.1230 
.1230 
.1230 
.1230 
.1230 
.1230 
. 1873 
.1873 
.1873 
.1873 
.1873 
.1873 
.1873 
.1873 
.1873 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.2300 
.3123 
.3123 
. 3123 
.3123 
. 3123 
. 3123 
.3123 
.3123 
. 3123 
. 3730 
.5730 
.3730 
.3730 
.3750 
.3750 
. 3750 
.3750 
. 3730 
.3730 
.3750 
.3750 
.3750 
3750 
3750 
3750 
3750 
4375 
4373 
4575 


• UUOO 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
. 0000 
.0000 
.0000 
.0000 
. 0000 
.0000 
. 0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
. 0000 
.0000 
.0000 
.0000 
,0000 
,0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 


. UUOO 
.0000 
.0000 
.0000 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 


. uuou 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
. 0000 
. 0000 
. 0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
,0000 
,0000 
oooo 
oooo 
0000 
oooo 
oooo 
oooo 
oooo 


99 

233 

2 

. 3750 

.4575 

.0000 

.0000 

1 00 

2 35 

2 

. 5000 

.4575 

.0000 

.0000 

101 

237 

2 

.6250 

.4375 

.0000 

.0000 

102 

239 

2 

.7500 

. 4375 

.0000 

.0000 



IU > 

24 1 

2 

104 

245 

2 

105 

245 

5 

106 

248 

2 

107 

250 

5 

108 

255 

2 

109 

255 

5 

110 

258 

2 

111 

260 

5 

112 

265 

2 

115 

265 

5 

114 

268 

2 

115 

270 

5 

116 

275 

2 

117 

275 

5 

118 

278 

2 

119 

280 

3 

120 

285 

2 

121 

285 

5 


IZZ 288 2 

125 290 2 

124 292 2 

125 294 2 

126 296 2 


127 

298 

2 

128 

500 

2 

129 

502 

2 

150 

504 

2 

151 

506 

3 

152 

509 

2 

155 

511 

3 

1 54 

514 

2 

155 

516 

3 

156 

319 

2 

157 

321 

3 

158 

324 

2 

159 

326 

3 

140 

329 

2 

141 

331 

3 

142 

334 

2 

145 

336 

3 

144 

339 

2 

145 

341 

3 

146 

344 

2 

147 

346 

3 

148 

349 

2 

149 

351 

2 

150 

353 

2 

151 

355 

2 

152 

357 

2 

155 

359 

2 

154 

361 

2 

155 

363 

2 

156 

365 

2 

157 

367 

3 

158 

370 

2 

159 

372 

3 

160 

375 

2 

161 

377 

3 

162 

380 

2 

165 

382 

3 

164 

385 

2 


.«/>U . 4 >/> 

1.0000 .4575 

.0000 .5000 

.0625 .5000 

•1250 .5000 

■1875 .5000 

.2500 .5000 

•3125 .5000 

*3750 .5000 

*^375 .5000 


.5000 

.5000 

.5625 

.5000 

.6250 

.5000 

.6875 

. 5000 

.7500 

.5000 

.8125 

. 5000 

.8750 

.5000 

.9375 

.5000 

1 .0000 

.5000 

.0000 

.5625 

.1250 

. 5625 

.2500 

.5625 

. 3750 

. 5625 

.5000 

.5625 

. 6250 

.5625 

.7500 

.5625 

.8750 

.5625 

1.0000 

.5625 

.0000 

.6250 

.0625 

.6250 

.1250 

.6250 

.1875 

.6250 

.2500 

.6250 

.3125 

.6250 

.3750 

.6250 

.4375 

.6250 

.5000 

.6250 

.5625 

.6250 

.6250 

.6250 

.6875 

.6250 

.7500 

. 6250 

.8125 

.6250 

.8750 

.6250 

.9575 

.6250 

l .0000 

.6250 

.0000 

.6875 

.1250 

.6875 

.2500 

.6875 

.3750 

.6875 

.5000 

.6875 

.6250 

.6875 

.7500 

.6875 

.8750 

.6875 

.0000 

.6875 

.0000 

.7500 

.0625 

.7500 

.1250 

.7500 

.1875 

.7500 

. 2500 

.7500 

.3125 

.7500 

. 3750 

.7500 

.4375 

.7500 


.uuuu 


. 0000 


.0000 

.0000 

.0000 


.0000 

.0000 

. 0000 


.0000 

.0000 

.0000 


.0000 

.0000 

.0000 


.0000 

.0000 

.0000 


.0000 

.0000 

.0000 


.0000 

.0000 

. 0000 


.0000 

.0000 

.0000 


.0000 

.0000 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

0000 

0000 

.0000 

0000 

0000 

.0000 

0000 

0000 

.0000 

0000 

0000 

.0000 


.uuuu 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

0000 

0000 

0000 

0000 

0000 

0000 


165 

387 

3 

. 5000 

.7500 

.0000 

.0000 

. 0000 

166 

590 

2 

. 5625 

.7500 

.0000 


. 0000 

167 

392 

3 

.6250 

.7500 

.0000 

.0000 

. 0000 

1 68 

395 

2 

. 6875 

.7500 

.0000 


. OQOO 



ib 9 

39 / 

> 

. />uu 

170 

400 

2 

. 8 1 2 > 

171 

402 

3 

. 87>0 

172 

403 

2 

. 937 > 

173 

407 

3 

1.0000 

174 

410 

2 

.0000 

173 

412 

2 

. 1 2>0 

176 

414 

2 

. 2>00 

177 

416 

2 

. 37>0 

178 

418 

2 

.>000 

179 

420 

2 

. 62>0 

180 

422 

2 

. 7 > 0 Q 

181 

424 

2 

. 87>0 

182 

426 

2 

1.0000 

183 

428 

3 

.0000 

184 

431 

2 

.0623 

183 

433 

3 

.1230 

186 

436 

2 

.1873 

187 

438 

3 

.2300 

188 

441 

2 

. 3123 

189 

443 

3 

.3730 

190 

446 

2 

.4373 

191 

448 

3 

.3000 

192 

431 

2 

.3623 

193 

433 

3 

.6230 

194 

436 

2 

.6873 

193 

438 

3 

.7300 

196 

461 

2 

.8123 

197 

463 

3 

.8750 

198 

466 

2 

.9375 

199 

468 

3 

1 .0000 

200 

471 

2 

.0000 

201 

473 

2 

.1250 

202 

473 

2 

.2500 

203 

477 

2 

.3750 

204 

479 

2 

.5000 

203 

481 

2 

.6250 

206 

483 

2 

.7500 

207 

483 

2 

.8750 

208 

487 

2 

1.0000 

209 

489 

3 

.0000 

210 

492 

2 

.0625 

211 

494 

3 

.1250 

212 

497 

2 

.1875 

213 

499 

3 

.2300 

214 

302 

2 

.3125 

213 

304 

3 

.3750 

216 

>07 

2 

.4375 

217 

>09 

3 

.>000 

218 

>12 

2 

.>625 

219 

>14 

3 

,6250 

220 

>17 

2 

.6875 

221 

>19 

3 

.7500 

222 

>22 

2 

.8125 

223 

>24 

3 

.8750 

224 

>27 

2 

.9375 

223 

>29 

3 

1.0000 


. />uu 

. uuuu 

. uuuu 

.7500 

.0000 


.7500 

.0000 

.0000 

.7500 

.0000 


.7500 

.0000 

.0000 

.8125 

.0000 


.8125 

.0000 


.8125 

.0000 


.8125 

.0000 


.8125 

. 0000 


.8125 

.0000 


.8125 

.0000 


.8125 

.0000 


.8125 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.8750 

.0000 


.8750 

.0000 

.0000 

.9375 

.0000 


.9375 

.0000 


.9375 

.0000 


.9375 

.0000 


.9375 

.0000 


.9375 

.0000 


.9375 

.0000 


.9375 

.0000 


.9375 

.0000 


1.0000 

.0000 

.0000 

1.0000 

.0000 


1.0000 

.0000 

.0000 

1.0000 

.0000 


1.0000 

.0000 

.0000 

1.0000 

.0000 


1 .0000 

.0000 

.0000 

1 .0000 

.0000 


1 .0000 

.0000 

.0000 

1 .0000 

.0000 


1.0000 

.0000 

.0000 

1.0000 

.0000 


1.0000 

.0000 

.0000 

1.0000 

.0000 


1.0000 

.0000 

.0000 

1.0000 

.0000 


1.0000 

.0000 

.0000 


. uuuu 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 
.0000 


BOUNDARY CONDITIONS 
******************* 


NODE 

U-FIXED PRESSURE 

V-F I XED 

209 


4.0000 

209 

5.0000 


210 


4.0000 



Cl u 

7 . UUUU 


211 


4.0000 

211 

5.0000 

212 


4.0000 

212 

5.0000 

21 ? 


4.0000 

215 

5.0000 

214 


4.0000 

214 

5.0000 

215 


4.0000 

215 

5.0000 

216 


4.0000 

216 

5.0000 

217 


4.0000 

217 

5.0000 


218 


4.0000 

218 

5.0000 


219 


4.0000 

219 

5.0000 


220 


4.0000 

220 

5.0000 


221 


4.0000 

221 

5.0000 


222 


4.0000 

222 

5.0000 


225 


4.0000 

223 

5.0000 


224 


4.0000 

224 

5.0000 


225 


4.0000 

225 

5.0000 


1 


.0000 

2 


.0000 

3 


.0000 

4 


.0000 

5 


.0000 

6 


.0000 

7 


.0000 

8 


.0000 

9 


.0000 

10 


.0000 

11 


.0000 

12 


.0000 

13 


.0000 

14 


.0000 

15 


.0000 

16 


.0000 

17 


.0000 

1 

2.0000 


18 

2.0000 


27 

2.0000 


44 

2.0000 


53 

2.0000 


70 

2.0000 


79 

2.0000 


96 

2.0000 


105 

2.0000 


122 

2.0000 


131 

2.0000 


148 

2.0000 


157 

2.0000 


174 

2.0000 



183 

2.0000 

200 

2.0000 

209 

2.0000 

17 

- 2.0000 



4b -4.U0UU 
43 - 2.0000 

52 - 2.0000 
69 - 2.0000 

78 - 2.0000 

95 - 2.0000 

104 - 2.0000 

121 - 2.0000 
150 - 2.0000 

147 - 2.0000 

156 - 2.0000 

173 - 2.0000 

182 - 2.0000 
199 - 2.0000 

208 - 2.0000 
225 - 2.0000 

1 .0000 



VALT VALUES AT 

GAUSS-POINTS 


1 

1.0 

1.0 

1.0 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

2 

1.0 

1.0 

1.0 

1 .0 

1.0 

1.0 

1.0 

1.0 

1 .0 

3 

1.0 

1.0 

1.0 

4.0 

1.0 

1.0 

4.0 

1.0 

1.0 

4 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

5 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1 .0 

6 

1.0 

1.0 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

1 .0 

7 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1 .0 

8 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1 .0 

1.0 

1.0 

9 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1 .0 

1.0 

1 .0 


10 4.0 4.0 4.0 4.0 1.0 1.0 1.0 1.0 1.0 

11 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

12 4.0 1.0 1.0 4.0 4.0 1.0 4.0 4.0 1.0 

13 1.0 1.0 1.0 4.0 1.0 1.0 4.0 1.0 1.0 

14 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 

15 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 

16 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 

17 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 

18 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1.0 1.0 

19 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

20 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

21 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

22 4.0 1.0 1.0 4.0 4.0 1.0 4.0 4.0 1.0 

23 1.0 1.0 1.0 4.0 1.0 1.0 4.0 1.0 1.0 

24 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 

25 4.0 4.0 4.0 4.0 1.0 1.0 1.0 1.0 1.0 

26 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

27 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

28 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

29 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

30 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

31 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

32 1.0 1.0 1.0 4.0 4.0 1.0 4.0 4.0 1.0 

33 1.0 4.0 4.0 1.0 4.0 4.0 1.0 1.0 1.0 

34 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

35 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

36 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

37 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

38 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

39 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

40 1.0 1.0 1.0 1.0 1.0 4.0 4.0 4.0 4.0 

41 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 

42 1.0 1.0 4.0 1.0 1.0 4.0 1.0 1.0 1.0 

43 1.0 4.0 4.0 1.0 4.0 4.0 1.0 1.0 4.0 

44 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 


45 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

46 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

47 1.0 1.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 

48 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 



49 

1 .0 

1 .0 

1 .0 

1 .0 

1.0 

1 .u 

i .U 

1.0 

1 . u 

50 

1.0 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

51 

1 .0 

1 .0 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

52 

1.0 

1.0 

4.0 

1.0 

1.0 

4.0 

1.0 

1.0 

1.0 

53 

1 .0 

4.0 

4.0 

1.0 

4.0 

4.0 

1,0 

1 .0 

4.0 

54 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

55 

1 .0 

1.0 

1.0 

1.0 

1.0 

4.0 

4.0 

4.0 

4.0 

56 

1 .0 

1 .0 

1.0 

1.0 

1.0 

1 .0 

1.0 

1.0 

1.0 

57 

1 .0 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

58 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

59 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

60 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

1 .0 

1.0 

1.0 

61 

1 .0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

62 

1.0 

1.0 

4.0 

1.0 

1.0 

4.0 

1.0 

1.0 

1.0 

63 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

64 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 


RESULTS 

FOR ITERATION 

NUMBER 

1 

NODE 

U-VELOCITY 1 

PRESSURE 

V-VELOCITY 

1 

. 2 00E + 01 

. 000E+00 

. 000E+00 

2 

. 1 74E+01 


. O0GE+00 

3 

.147E+01 - 

. 176E+03 

. QQ0E+00 

4 

.12 1E+01 


. 000E+00 

5 

. 695E + 00 - 

. 1 29E+04 

. 000E + Q0 

6 

. 21 7E+0G 


. OOOE+OO 

7 

- . 245E+00 - 

. 118E+04 

. 000E+00 

8 

- . 671 E+00 


.OOOE+OO 

9 

- . 972E+00 

. 241E+03 

.000E+00 

10 

123E+01 


.OOOE+OO 

11 

- . 1 45E+G1 

. 596E+02 

.OOOE+OO 

12 

- . 1 64E+01 


.OOOE+OO 

13 

178E+01 

. 284E+0 3 

.OOOE+OO 

14 

- . 1 87E+01 


.OOOE+OO 

15 

- . 1 94E+01 

. 484E+0 3 

.OOOE+OO 

16 

- . 197E+01 


.OOOE+OO 

17 

- . 200E + 0 1 

. 565E+03 

.OOOE+OO 

18 

. 200E + 01 


. 281 E + 00 

19 

. 1 46E+0 1 


. 262E + 00 

20 

. 0Q0E+00 


.OOOE+OO 

21 

- . 1 68E+00 


. 41 8E+00 

22 

- . 925E+00 


. 294E+00 

23 

- . 143E+01 


. 21 4E+00 

24 

- . 176E+01 


. 1 1 5E+00 

25 

- . 193E+01 


. 490E-01 

26 

- . 200E+0 1 


. 267E-Q1 

27 

. 200E+0 1 

. 250E + Q2 

. 574E+00 

28 

. 165E+01 


. 578E+00 

29 

. 1 36E+0 1 

. 148E+03 

. 630E+00 

30 

. 867E+0Q 


. 784E + 00 

31 

. 450E+00 - 

. 531E+G3 

. 872E+00 

32 

.223E+00 


. 861E+00 

33 

- . 294E-0 1 - 

. 589E+0 3 

. 835E+00 

34 

- . 423E+00 


. 707E+00 

35 

-.791E+00 - 

.439E+Q3 

. 552E+00 

36 

- . 1 08E+01 


. 493E+00 

57 

- . 1 34E+01 - 

. 201E+02 

. 407E+00 

38 

- . 1 55E+01 


. 326E+00 

39 

- . 171E+01 

. 423E+03 

.259E+00 

40 

-.182E+01 


. 1 80E+00 

41 

- . 1 91E+01 

. 571E+03 

. 1 1 6E+00 


42 

- . 1 96E+01 


. 704E-01 

43 

- . 200E + 0 1 

. 605E+03 

. 526E-G1 

44 

. 200E + 01 


. 906E + 00 

45 

. 1 2 1 E + 01 


. 1 07E+01 



46 

.4>6t+UU 

47 

- . 881E-Q1 

48 

- . 605E+00 

49 

-.116E+Q1 

*0 

- . 1 62E+01 


-.188E+01 

52 

- . 200E+01 

53 

. 200E+01 

54 

. 1 58E+01 

55 

. 111E+01 

56 

. 726E+00 

57 

. 456E+00 

58 

. 1 32E+00 

59 

- . 1 54E+00 

60 

- . 368E+00 

61 

** . 608E+00 

62 

- . 806E+0Q 

63 

- - 990E+00 

64 

-.127E+01 

65 

-.151E+01 

66 

- . 170E+01 

67 

- . 1 83E + 01 

68 

- . 1 94E+01 

69 

- . 200E + 01 

70 

. 200E+01 

71 

. 1 07E + 01 

72 

. 449E+00 

73 

- . 1 44E+00 

74 

- . 625E+00 

75 

- . 1 03E+01 

76 

- . 140E+01 

77 

- . 175E+01 

78 

- . 200E+0 1 

79 

. 2GGE+Q1 

80 

. 153E+01 

81 

. 1 07E+01 

82 

. 779E+00 

83 

. 458E+00 

84 

. 147E+00 

85 

- . 1 34E+00 

86 

- . 395E+00 

87 

- . 639E + 00 

88 

- . 868E + 00 

89 

- . 1 07E+01 

90 

- . 1 25E+01 

91 

-.143E+01 

92 

- . 1 57E+Q1 

93 

-.170E+01 

94 

“ . 1 87E+01 

95 

- . 200E+01 

96 

. 200E+01 

97 

. 1 14E+01 

98 

.481E+00 

99 

-.116E+00 

100 

- . 6 35E+00 

101 

-.109E+01 

102 

- . 146E+0 1 

103 

- . 176E+01 

104 

- . 200E+01 

105 

. 200E+01 

106 

. 1 51 E+0 1 

107 

. 1 1 9E+01 


108 

. 809E + 00 

109 

. 482E+00 

110 

. 1 90E+00 

111 

- . 905E-01 


-427E+03 

. 1 1 /t + Ui 
.106E+Q1 
. 943E+00 
. 672E+00 
. 405E+00 
. 192E+00 
. 896E-01 
. 1 35E+01 

. 645E+0 3 

. 1 37E+01 
. 1 47E+01 

. 770E + 0 3 

. 1 52E+01 
.145E+01 

. 339E+03 

. 1 38E+01 
. 1 30E+01 

640E+03 

. 1 24E+01 
.116E+01 

412E+03 

. 1 08E+01 
. 1 OOE+Ol 

382E+03 

. 782E+00 
. 539E+00 

658E+03 

. 398E+Q0 
. 270E+00 

878E+03 

. 1 81E+00 
. 1 49E+Q0 

111E+04 

. 1 7 5E+01 
. 1 93E+0 1 
. 1 75E + 01 
.158E+01 
. 1 38E+01 
.114E+01 
. 82 3E+00 
. 426E+00 
. 223E+00 
. 2 1 4E + 01 

728E+02 

. 229E+01 
. 2 33E+01 

360E+03 

. 2 1 9E+01 
, 206E + 01 

114E+03 

. 1 9 5E + 01 
. 1 8 5E + 01 

502E+03 

. 1 74E+01 
.161E+01 

125E+04 

. 1 47E + 01 
.131E+01 

151E+04 

. 1 1 5E+01 
. 973E+00 

118E+04 

. 807E+00 
. 657E + 00 

120E+04 

. 407E + 00 
. 282E + 00 

232E+04 

. 273E + 01 
. 26 5E + 01 
. 2 37E + 01 
. 21 2E+01 
. 1 85E+01 
. 1 54E+01 
. 1 1 4E+01 
. 726E+00 
. 440E + 00 
. 329E+01 

135E+03 

. 320E+01 
. 304E+01 

557E+00 

. 287E + 01 
. 269E + 01 

671 E+0 1 

.254E+01 
. 2 39E + 01 



112 

- . 555L+UU 

115 

-.603E+00 

114 

- , 845E + 00 

115 

- . 1 08E+01 

11 6 

- . 1 30E+01 

117 

-.150E+01 

118 

- . 1 68E + 01 

119 

- . 1 83E+01 

120 

- . 1 90E+Q1 

121 

- . 200E + 01 

122 

. 200E+01 

123 

.111E+01 

124 

. 465E+00 

125 

- . 567E-01 

126 

- . 548E+00 

127 

- . 1 03E+01 

128 

- . 1 51E+01 

129 

- . 1 86E+01 

130 

- . 2Q0E+01 

131 

. 200E+0 1 

132 

, 1 50E+Q1 

133 

. 997E+Q0 

134 

. 694E+00 

135 

. 442E+00 

136 

. 1 99E+00 

137 

- . 1 07E-0 1 

138 

- . 224E+00 

139 

- . 456E+00 

140 

- . 691E+00 

141 

- . 939E+00 

142 

- . 1 21E+01 

143 

-.147E+01 

144 

- . 171E+0 1 

145 

-.188E+01 

146 

- . 205E+01 

147 

- . 20QE+0 1 

148 

. 200E+0 1 

149 

. 1 1 1E+01 

150 

. 459E+00 

151 

. 449E-0 1 

152 

- . 305E + 00 

153 

- . 775E+00 

154 

- . 1 39E+0 1 

155 

“ . 1 90E+01 

156 

“ . 200E+01 

157 

. 200E+0 1 

158 

. 1 74E+0 1 

159 

. 1 30E + 01 

160 

. 1 01E+01 

161 

. 814E+00 

162 

. 484E+00 

163 

. 1 74E+00 

164 

. 447E-0 1 

165 

- . 984E-01 

166 

- . 286E+00 

167 

- . 502E+00 

168 

-.781E+00 

169 

-.118E+01 

170 

- . 1 62E+01 

171 

- . 1 94E+01 

172 

-.178E+01 

173 

- . 200E+0 1 


174 

. 200E + 0 1 

175 

. 1 67E+01 

176 

.136E+01 

177 

. 979E+00 


. 472E+03 

- 22>t+Ui 
-210E+01 

. 907E+03 

. 1 9 3E+01 
. 1 7 5E+01 

. 1 62E+04 

.155E+01 
. 1 32E+01 

. 2 3 5E + 04 

. 1 07E+01 
. 832E+00 

. 299E+04 

. 655E+0Q 
. 500E+00 

- . 1 31E+04 

. 376E+01 
. 336E+01 
. 300E+01 
. 265E+01 
. 2 34E+01 
.198E+01 
. 1 52E+01 
. 992E+00 
. 549E+O0 
. 435E+01 

- . 1 36E + 04 

. 403E+01 
. 362E+01 

-.466E+03 

. 343E+01 
. 325E+01 

~ . 438E+02 

. 306E+01 
. 289E+01 

. 1 45E+03 

. 273E+01 
. 258E+01 

. 878E+03 

. 242E+01 
. 224E+01 

. 1 65E+04 

. 201E+01 
. 1 75E+0 1 

. 289E+04 

. 1 45E+01 
.112E+01 

•455E+04 

. 781 E+00 
. 566E+00 

- . 252E+04 

. 473E+01 
.419E+01 
. 346E+01 
. 308E+01 
. 280E+G1 
-250E+01 
. 203E+0 1 
. 126E+01 
. 597E+00 
. 507E+0 1 

-.231E+04 

. 494E+0 1 
. 461 E+Q 1 

- . 1 86E+04 

. 421E+01 
. 386E+01 

- . 850E+03 

. 349E+01 
. 321E+01 

- . 2 96E+03 

. 309E+01 
. 298E+01 

. 245E+03 

. 287E+0 1 
. 276E+01 

. 1 61 E+04 

. 262E+0 1 
. 2 37E + 01 

. 527E+04 

. 206E+01 
. 1 52E+0 1 

. 5 1 5E+04 

. 863E+00 
. 61 7E+00 


.537E+01 
. 478E+01 
. 398E+01 
. 355E+01 


178 

. 257E+00 


. 314E+01 

179 

-.11 7E+00 


. 300E+01 

180 

-.751E+00 


. 282E+01 

181 

- . 1 84E+01 


. 1 78E+01 

182 

- . 200E+01 


. 535E+00 

183 

. 200E+01 

- . 31 6E+04 

. 60 3E+G1 

184 

. 1 53E+Q1 


.541E+01 

185 

. 226E+01 

- . 373E+04 

. 440E+01 

186 

. 242E+01 


. 408E+01 

187 

. 223E+01 

- . 1 31E+04 

. 416E+01 

188 

. 223E+01 


. 405E+01 

189 

. 21 5E+01 

150E + 04 

. 379E+01 

190 

. 1 83E+01 


. 358E+01 

191 

. 1 56E+01 

- . 754E + 0 3 

. 354E+01 

192 

. 1 1 3E+01 


. 349E+01 

193 

. 51 6E+00 

- . 872E+02 

. 332E+01 

194 

. 21 0E+00 


. 31 1 E+01 

195 

- . 832E-01 

- . 748E+Q 3 

. 306E+01 

196 

- . 300E+00 


. 347E+01 

197 

- . 108E+01 

.679E+04 

. 309E+01 

198 

- . 300E+0 1 


. 7 1 3E+00 

199 

- . 200E+0 1 

. 55 3E+04 

- . 822E+00 

200 

. 200E+01 


. 621 E+01 

201 

. 267E+0 1 


. 384E+01 

202 

. 382E+0 1 


. 425E+01 

203 

. 327E+0 1 


. 390E+01 

204 

. 345E+0 1 


. 380E+01 

205 

. 207E+0 1 


. 392E+01 

206 

. 1 94E+01 


. 322E+01 

207 

- . 41 5E+00 


. 463E+01 

208 

“ . 200E+0 1 


127E+01 

209 

. 200E+0 1 

- . 1 83E+0 5 

. 400E+Q1 

210 

. 500E+Q 1 


. 400E+01 

211 

. 500E+01 

. 386E+04 

. 400E+01 

212 

. 500E+0 1 


. 400E+01 

213 

. 500E+01 

- . 461 E+04 

. 4Q0E+01 

214 

. 500E+0 1 


. 400E+01 

215 

. 500E+-01 

- . 1 01E+04 

. 400E+01 

216 

. 500E+0 1 


. 400E+0 1 

217 

. 500E+01 

- . 498E+01 

. 400E+01 

218 

. 500E+01 


. 400E+01 

219 

. 500E+01 

- . 51 9E+04 

. 400E+01 

220 

. 500E+01 


. 400E+01 

221 

. 500E+01 

. 669E+Q4 

. 400E+01 

222 

. 500E4-01 


. 400E+01 

223 

. 500E+0 1 

- . 1 38E+05 

. 400E+01 

224 

. 500E+01 


. 400E+01 

225 

- . 2Q0E+01 

. 42 9E+05 

. 400E+01 


LARGEST CHANGE OCCURS AT DOF NUMBER 531 CHANGE = 


1.0000 



RESULTS FOR ITERATION NUMBER 24 


NODE 

U-VELOC I TV 

PRESSURE 

V-VELOCITY 

1 

. 2Q0E+0 1 

. OOOE+OO 

. OOOE+OO 

2 

.187E+01 


. OOOE+OQ 

3 

. 1 76E+0 1 - 

. 265E + 02 

.OOOE+OO 

4 

. 166E+01 


.OOOE+OO 

5 

. 1 08E+0 1 - 

. 144E+03 

.OOOE+OO 

6 

. 691E+00 


.OOOE+OO 

7 

. 264E+00 - 

.193E+03 

.OOOE+OO 

8 

- . 283E+0Q 


.OOOE+OO 


9 

- . 600E + 00 

- . 1 1 4E + 0 3 

.OOOE+OO 

10 

- . 824E + 00 


.OOOE+OO 

11 

-.103E+01 

- . 889E + 0 2 

.OOOE+OO 

12 

- . 1 28E + 01 


.OOOE+OO 

1? 

- . 146E+01 

- . 752E + 02 

.OOOE+OO 

14 

-.164E+01 


•OOOE+OO 

15 

-. 179E+01 

- . 732E+02 

.OOOE+OO 

16 

- . 189E + 01 


•OOOE+OO 

17 

- . 200E+01 

- . 844E + 02 

•OOOE+OO 

18 

. 200E+0 1 


. 143E+00 

19 

.OOOE+OO 


•OOOE+OO 

20 

. 943E + 00 


. 464E+00 

21 

. 278E + 00 


•442E+00 

22 

- . 564E+00 


. 297E+00 

23 

102E + 01 


•235E+00 

24 

145E + 01 


. 170E+00 

25 

- • 179E + 01 


•120E+00 

26 

- . 200E+01 


. 110E+00 

27 

. 200E+01 

- . 1 57E+02 

. 287E+00 

28 

, 177E + 01 


. 306E+00 

29 

. 1 62E + 01 

. 193E + 02 

. 380E+00 

30 

. 1 08E+0 1 


. 710E+00 

31 

. 677E+00 

- . 740E+02 

. 880E+00 

32 

. 508E+00 


•881E+00 

35 

. 329E+00 

- . 1 34E + 03 

. 867E+00 

34 

- . 395E-01 


. 699E+00 

55 

-.451E+00 

- . 1 6 1 E + 0 3 

. 487E+00 

36 

- . 702E+00 


. 468E+00 

37 

- . 976E+00 

- . 1 20E + 03 

. 447E+00 

38 

- . 1 24E+01 


•425E+00 

39 

-. 142E+01 

- . 575E + 02 

•413E+00 

40 

- . 159E + 01 


. 349E+00 

41 

- . 176E+01 

- . 567E + 02 

. 278E+00 

42 

- . 1 88E + 01 


. 236E+00 

43 

- . 200E+0 1 

-.771E+02 

•219E+00 

44 

. 200E + 01 


. 51 1E + 00 

45 

. 1 34E + 01 


•826E+00 

46 

. 633E+00 


. 1 07E+01 

47 

. 229E + 00 


. 1 01E+0 1 

48 

- . 256E+00 


. 947E+0 0 

49 

- . 799E+00 


•743E+00 

50 

- . 1 37E+0 1 


. 563E+00 

51 

- . 173E+01 


. 416E+00 

52 

“ . 200E+0 1 


. 341E+00 

53 

. 200E+0 1 

. 404E+02 

.915E+00 



54 

. 1 62E+01 


. 954E + 00 

55 

. J17E+01 

. 732E + 02 

. 121E+01 

56 

. 829E+00 


•133E+01 

57 

. 622E+00 - 

. 930E+02 

.130E+01 

58 

. 358E+00 


. 1 25E+01 

59 

. 1 38E+00 

. 264E + 02 

. 121E + 01 

60 

- . 666E-01 


-118E+01 

61 

- . 289E+00 

. 21 8E + 02 

. 1 15E+01 

62 

- . 477E+00 


.114E+01 

65 

-.637E+00 - 

. 605E+02 

.112E+01 

64 

- . 950E+00 


. 926E+00 

65 

- . 1 25E+01 - 

. 937E+02 

. 697E+O0 

66 

- . 148E+01 


. 628E+00 

67 

- . 169E+01 

. 639E+02 

. 545E+00 

68 

- . 1 87E+01 


. 480E+00 

69 

- . 200E+01 

. 345E+02 

. 472E+00 

70 

.2Q0E+01 


. 125E+01 

71 

. 1 G9E+Q1 


. 1 64E+01 

72 

. 590E+00 


. 1 54E+Q 1 

73 

. 821E-01 


. 146E+01 

74 

“ . 343E+00 


.134E+01 


75 

- . 714E + 0Q 


. 1 22E+0 1 

76 

-.111E+01 


.109E+01 

77 

- . 1 58E + G1 


. 793E+00 

78 

- . 200E+01 


. 605E+Q0 

79 

. 200E + 01 

. 1 33E+03 

.160E+01 

80 

. 1 52E+01 


. 1 86E+0 1 

81 

• 1 06E + 0 1 

. 873E+01 

-198E+01 

82 

. 8 38E+00 


.188E+01 

83 

. 566E+00 

- . 559E+02 

. 181E+01 

84 

. 292E+00 


. 1 75E+0 1 

85 

. 524E-01 

“ . 298E+02 

. 1 69E+01 

86 

- . 176E+00 


.163E+01 

87 

- . 396E+00 

- . 320E+01 

. 1 55E+01 

88 

- . 603E + 00 


.147E+01 

89 

-.793E+00 

. 985E+02 

. 1 39E+01 

90 

- . 967E+00 


. 1 32E+01 

91 

-.115E+01 

. 1 02E+03 

. 1 25E+01 

92 

-.131E+Q1 


. 1 20E + 0 1 

93 

- . 147E+01 

- . 878E+01 

.115E+01 

94 

-.177E+01 


. 923E + 00 

95 

- . 200E + 01 

.141E+02 

. 745E + 00 

96 

. 200E + 01 


. 228E + 01 

97 

. 1 12E+01 


. 227E+01 

98 

. 550E+00 


. 208E+01 

99 

. 278E-01 


. 1 92E+0 1 

100 

- . 429E+00 


.177E+01 

101 

- . 839E+00 


.161E+01 

102 

~ . 120E+01 


. 1 42E+01 

103 

- . 1 55E+01 


. 1 25E+01 

104 

- . 200E+01 


. 1 06E + 01 

105 

. 200E + 01 

. 265E+03 

. 295E+01 

106 

. 1 45E + 0 1 


. 280E+01 

107 

.116E+01 

. 357E+02 

. 265E+01 

108 

.800E+00 


. 25GE+01 

109 

. 509E+00 

- . 1 74E+02 

. 235E + 01 

110 

. 254E+00 


. 224E+01 

111 

. 1 08E-0 1 

- . 592E+02 

. 21 5E+01 

112 

- . 21 8E+00 


. 207E+01 

113 

- . 436E+00 

-.11 6E+02 

. 1 99E+01 

114 

- . 651E+00 


. 1 90E+01 

115 

- . 857E + 00 

. 282E+02 

. 181E+01 

116 

- . 1 06E+01 


. 1 70E+01 

117 

-.125E+01 

. 957E+02 

.160E+01 

118 

-.143E+01 


. 1 50E+01 

119 

- .162E+01 

.168E+03 

. 140E+01 


10+301** 

1 0+3ZZZ * 

IO+300Z* 

18 r 

10+399** 


1 0+3/9 T * 

*8r 

10+36/** 

1 0+39Z1 * - 

r o+3ooz * 

181 

TO+390Z* 


I0+300Z * - 

zs r 

T0+3Z9Z* 


ro+399r *- 

ra r 

IQ+39ZZ * 


00+3168*- 

081 

I0+3t8Z* 


00+3911 *- 

6 /r 

I0+388Z ' 


00 + 31/1 *- 

8 /r 

ro+azu * 


00+31 Z9 * 

zzi 

10+3861 * 


ro+39*r* 

9/r 

I0+3Z1** 


ro+aosr * 

i/r 

10+309** 


IO+300Z* 

*zr 

IQ+310Z * 

10+3111* 

I 0+300Z * - 

izr 

T0+360Z * 


ro+3i/r *- 

ZZ T 

10+30ZZ * 

10+3011* 

ro+3r/r*- 

III 

I0+3*1Z* 


10+39*1*- 

dll 

10+38 trZ 

Z0+3T6Z*- 

ro+3zir *- 

691 

TO+309Z* 


00+39 18 * - 

891 

10+31 9Z * 

10+391/ * 

00+3119*- 

Z9I 

10+369Z * 


00+31/**- 

991 

tO+31 ZZ ‘ 

Z0+3I8Z*- 

00+36Z1 * - 

!9r 

I0+3ZZZ* 


QQ+310Z * - 

*91 

T0+318Z* 

10+3111*- 

r0-3808*- 

191 

10+31 Z1 * 


00+3181* 

Z9 1 

ro+31/i* 

10+3ZIZ*- 

00+39Z6 ’ 

19 r 

ro+3*o** 


ro+3111* 

091 

T0+3ZZ*' 

1 0+31*Z * - 

ro+3*9i* 

611 

ro+az*** 


10 + 3161 * 

811 

TO+38*** 

10+3Z9T * - 

IO+300Z * 

Z1I 

rO+3/61* 


T0+300Z * - 

9i r 

10+3061* 


ro+3i9r*- 

111 

I0+31ZZ * 


T0+3ZZ I * - 

*11 

I0+3**Z* 


Q0 + 36ZZ * - 

111 

T0+361Z* 


00+3191*- 

zir 

I0+39ZZ * 


tO-3111 ' - 

m 

r0+3901* 


00+360** 

on 

TO+3181 * 


T0+39Zr * 

6* t 

10+301* * 


T0 + 300Z * 

8*1 

r0+308t* 

10+36 T1 * 

I0+30OZ * - 

/* t 

T0+39ZT ’ 


10 + 318 r * - 

9*r 

ro+3izr 

ZO+3916* 

10+3091*- 

i*t 

I0+318T * 


ro+3i*r- 

**r 

10+3861* 

10+3/01* 

I0+39Z I ’ - 

i*i 

ro+artz* 


10+3*01*- 

Z*I 

I0+31ZZ* 

ro+3o/** 

00+3608*- 

i*r 

I0+3Z1Z * 


00+3/61 * - 

0*1 

I0+3I*Z* 

Z0+36 t* * - 

00+3161*- 

611 ! j 

10+3I1Z * 


00+3661*- 

811 

I0+3I9Z * 

Z0+30IZ* 

10-31 IZ * - 

/II 

ro+3izz* 


00+3991* 

911 

10+318Z ‘ 

Z0+3101 * - 

00+3/61* 

lit 

10+3Z6Z* 


00+31*9* 

*11 

10+3*11* 

10+39ZT *- 

00+3Z1 6 ’ 

111 

I0+3Z91 * 


r0+31!T* 

Z1I 

10+310** 

10+38H *- 

TO + 300Z * 

lit 

ro+3zir* 


ro+aooz*- 

Oil 

I0+3r9I* 


t0+3Z9I*- 

6ZI 

10+38/ t* 


T0+3/ZI *- 

8Z I 

lo+aroz* 


00+3118*- 
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